perm filename JAMLIB.JAM[UP,DOC]3 blob
sn#476181 filedate 1979-09-24 generic text, type C, neo UTF8
COMMENT ⊗ VALID 00016 PAGES
C REC PAGE DESCRIPTION
C00001 00001
C00003 00002 USING JAMLIB
C00005 00003 LIST OF DECLARATIONS (without comment)
C00013 00004 NUMERICAL ANALYSIS
C00014 00005 NUMERICAL ANALYSIS - SIMULTANEOUS EQUATIONS
C00027 00006 NUMERICAL ANALYSIS - CURVE FITTING
C00031 00007 NUMERICAL ANALYSIS - OPTIMISATION
C00039 00008 NUMERICAL ANALYSIS - MISCELANEOUS
C00044 00009 SIGNAL PROCESSING
C00045 00010 SIGNAL PROCESSING - FAST FOURIER TRANSFORMS
C00059 00011 SIGNAL PROCESSING - FILTERING
C00080 00012 SIGNAL PROCESSING - LINEAR PREDICTION
C00100 00013 SIGNAL PROCESSING - SOUND FILE IO
C00122 00014 DISPLAY
C00172 00015 MISCELANEOUS - INREAL, ININT, DEFAULT
C00175 00016 MISCELANEOUS - LINE SEGMENT FUNCTION IO
C00193 ENDMK
C⊗;
USING JAMLIB
Binary library: DSK:JAMLIB.REL[SUB,SYS]
Sources: MUSIC0:[LIB,JAM]
Compile list: MUSIC0:JAMCM1.CCL[LIB,JAM], JAMCM2.CCL
Library list: MUSIC0:JAMRL1.CCL[LIB,JAM], JAMRL2.CCL
Complete List of Calls: DSK:JAMLIB.SAI[SUB,SYS], and page 3
of this document
JAMLIB is a library of SAIL or FAIL callable routines for doing numerical
analysis, sound processing, and display. The complete document list may
be found on JAMLIB.SAI[SUB,SYS]. This document gives only a brief
outline of the routines present and of how to read the documentation, as well
as tutorial examples in the use of some of the programs.
LIST OF DECLARATIONS (without comment)
Take a look at JAMLIB.SAI[SUB,SYS] and possibly JAMFOR.SAI[SUB,SYS] for help.
The sources used in compiling these libraries exist on the MUSIC0 User Disk Pack
under the PPN [LIB,JAM]. Below are the procedures
defined in JAMLIB without the documentation on SUB,SYS.
(from JAMLIB.HDR[SUB,SYS])
require "JAMLIB.REL[SUB,SYS]" library;
Comment - Numerical Analysis routines;
External real procedure EVALP(real array P; Integer N; real x);
External procedure LINFIT(real array F,Slope,Intercept,RMS;
integer N,width);
External integer procedure LSOLVE(real array A,H,F;
integer M,N,P);
External integer procedure DCOMP1(integer M,N; real array qr,alpha;
integer array pivot);
External integer procedure SOLVE1(integer M,N; real array qr,alpha;
integer array pivot; real array r,y);
External integer procedure LSQ(real array X,F,A; integer M,N);
External integer procedure RSQ(real array F,A; integer M,N; real SX,DX);
External Real procedure PMIN(real y0,y1,y2,y3; reference real minval);
External Real procedure PMAX(real y0,y1,y2,y3; reference real minval);
External integer procedure ROOTER(real array A; integer N;
Reference real array Xr,Xi; real tolerance; integer maxit);
External integer procedure DECOMP(integer N; real array A,LU);
External procedure SOLVE(integer N; real array B,X);
External integer procedure IMPROVE(integer N; real array A,B,X);
External integer procedure QSOLVE(integer N; real array A,X,B);
External integer procedure INVERT(integer N; real array A,INV);
Comment - Signal processing routines ;
External procedure INIALM(real array H,K,A; integer M; real fund,clock);
External procedure ALLBUM(real array X,Y,H; integer N);
External procedure ARG(real array REALPT,IMAG,MAGNITUDE,ANGLE;
integer N,lapnum; boolean pickup);
External real SAVSN,SAVCS,SAVAR;
External procedure ARTRAN(real array X,A,B; integer N; boolean inverse);
External procedure AUTOC(real array X; integer N;
real array R; integer M);
External procedure FFTARB(real array A,B; integer N,Nv,Ks);
FORTRAN PROCEDURE FRXFM;
External boolean procedure GFINI(real array H; real Wc1,Wc2,clock;
integer type,class,order);
External procedure GFILT(real array X,Y,H; integer N; boolean backwards);
External procedure GTRAN(real array X,A,B; Integer N,M,sN;
real phi,omega,clock; boolean pickup);
External real savla,savlb;
External integer savspt;
External procedure INVFLT(real array R,A,TEMP1,TEMP2; integer M);
External procedure RTRAN(real array X,A,B; integer log2N; boolean inverse);
External procedure SNDIN(integer channel,sample_number,N;
real array H,X; reference integer iflag);
External procedure SNDOUT(integer channel,N;
real array H,X; reference integer oflag);
External procedure WEIGHT(real array X; integer N;
real array D; integer M);
External procedure XRTRAN(real array A,B; integer N; boolean evaluate);
Comment - Display routines - ;
External procedure DSETUP(integer nwds; reference integer id);
External boolean procedure DRELS(reference integer id);
External procedure WRITE(integer id,pog);
External boolean procedure DWRITE(integer id,chan);
External procedure BUFCLR(integer id,nwds);
External procedure RVECT(integer id,dX,dY);
External procedure RIVECT(integer id,dX,dY);
External procedure AVECT(integer id,X,Y);
External procedure AIVECT(integer id,X,Y);
External procedure RPT(integer id,dX,dY);
External procedure APT(integer id,X,Y);
External procedure DTEXT(integer id; string text; real scale,angle);
External integer frame;
External procedure DAXIS(integer id; reference integer frame;
string xlabel,ylabel; reference real params);
External real params;
External procedure DAXISF(integer id; reference integer frame;
integer xlabel,ylabel; reference real params);
External procedure DPY(integer id; reference real params;
integer npts; reference real X);
External procedure DPYS(integer id; real array X;
integer N,Nslices,sliceN; real vxmin,vxmax);
External procedure DPYSL(integer id; real array X;
integer N,Nslices,sliceN; string xlabel,ylabel;
real vxmin,vxmax);
External string procedure ASC2S(integer text);
External procedure DPYSLF(integer id; real array x;
integer N,Nslices,sliceN; real vxmin,vxmax;
integer xlabel,ylabel);
External procedure AXIS(integer ID; real VMIN,VMAX;
reference real SCALE,OFFSET;
integer POS,MIN,MAX;
boolean XAXIS);
External procedure ASCALE(real VMIN,VMAX; reference real SCALE,OFF;
integer MIN,MAX;
reference integer TSTART,LTSTART,LTWID,BTWID;
reference real SI,DI);
External procedure XAXIS(integer ID,YPOS,XMIN,XMAX,
TSTART,LTSTART,LTWID,BTWID;
real SI,DI);
External procedure YAXIS(integer ID,XPOS,YMIN,YMAX,
TSTART,LTSTART,LTWID,BTWID;
real SI,DI);
External procedure TYPLOC(integer ymin,ymax);
NUMERICAL ANALYSIS
There are 3 branches of numerical analysis covered by the
routines in JAMLIB, plus a miscellaneous collection of useful odds and ends:
1) Simultaneous equations (matrix inversion)
2) Curve Fitting (splines, polynomials)
3) Optimisation (linear or non-linear)
4) Miscellaneous (Bessel functions, random numbers, etc)
We will try to discuss them one at a time.
NUMERICAL ANALYSIS - SIMULTANEOUS EQUATIONS
Given a system of simultaneous equations in the following form:
A[1,1]*X[1] + A[1,2]*X[2] + . . . + A[1,N]*X[N] = B[1]
A[2,1]*X[1] + A[2,2]*X[2] + . . . + A[2,N]*X[N] = B[2]
. . . . . . . . .
A[N,1]*X[1] + A[N,2]*X[2] + . . . + A[N,N]*X[N] = B[N]
where A[1:N,1:N] is a known NxN matrix and B[1:N] is a known vector
of length N, solve for X[1:N]. The routine that does this is called
QSOLVE and has the following form:
External integer procedure QSOLVE(integer N; real array A,X,B);
For instance, we can use this to do linear Gauss/Newton/Lagrange polynomial
fits in the following manner: If we have a data vector of observations
OBS taken at certain times TIMES, we might use the following program
to compute the coefficients of a polynomial that exactly matches the
points as follows:
for i←0 step 1 until N-1 do
for j←0 step 1 until N-1 do
A[i+1,j+1] ← TIMES[i+1]↑j;
err←qsolve(N,A,X,OBS);
The coefficients of the polynomial will now be in X. X[1] will be the
constant term, X[2] will be the coefficient of TIME↑1, X[3] of
TIME↑2, etc.
The QSOLVE routine returns an integer reporting the success or failure of
the computation. This integer should be interpreted as follows:
0 - Success. Correct solution vector in X.
1 - Matrix with zero row (singular)
2 - Singular matrix, about to divide by zero doing Gaussian elimination
3 - No convergence in IMPROVE, Matrix is nearly singular (eigenvalue
spread is too wide)
Error codes 4 through 7 mean that the routines could not obtain enough
scratch memory space to do the intermediate calculations and in short,
your matrix is too large.
There is the following routine for matrix inversion:
External integer procedure INVERT(integer N; real array A,INV);
You give it a square matrix A of dimensions A[1:N,1:N], and it returns
you the inverse of this matrix (if possible) in INV[1:N,1:N]. It also
returns the same error codes as QSOLVE (0=success, non-zero = failure).
Another possible usage for this kind routine is not for exact polynomial
fits but for linear least-squares fits. For this special case, there
is a more accurate way to do it.
The usual linear least-squares problem can be stated as follows:
Given that you have a M observations OBS[i] of some process at M
times that we will call TIMES[i], and you want to approximate these
observations with a weighted sum of N functions G[j] with weights H[j],
this can be done by minimizing the following quantity over the M
observations:
M N
E = ↔ (OBS[i] - ↔ H[j]*G[j](X[i]) )↑2
i=1 j=1
This is the sum of the squared differences between the observed values
at the M points and the approximate values from the calculation of
the weighted sum of the N functions at the given times of observation.
For example, if the functions are polynomials, then
(j-1)
G[j](X) = X
To solve this system for the unknown coefficients H[j], you do the following:
Make up a matrix A[1:M,1:N] where there are M data points (X[i]) and
N coefficients (H[j]) as follows:
A[i,j] ← G[j](X[i])
And make up the matrix of observations as taking the second index to
be always one. Then call the following routine:
External integer procedure LSOLVE(real array A,H,F; Integer M,N,P);
where M is the number of data points and N is the number of coefficients.
The routines actually solves the problem P times for P different
sets of observed values. Normally, one would set P=1 if one has only one
set of observations.
The routine returns an integer which will be zero upon success and one if the
matrix is singular. It returns 2 if there is not enough memory for
the internal calculations.
Notice that the formulation of several entirely different problems leads
to the same equations, such that this same routine can be used in many
different forms. This routine is good for any linear least-squares
problem. For instance, the functions G[j] might be functions of several
variables, not just one variable. It doesn't matter, because the X[i]
never appear explicitly, only implicitly as their function values,
G[j](X[i]) and the observed values. Thus, X[i] can be a vector, and the
functions G[j] can be functions of many variables. This says that there
are many different data vectors (rather than data points) that we wish to
minimize over.
Thus given the M observations in OBS[1:M] and the M times of the
observations in TIMES[1:M], the following SAIL program fragment will
calculate the coefficients H[j] of a J-1 order polynomial that
approximates the observations:
for i←1 step 1 until M do
for j←1 step 1 until N do
A[i,j] ← TIMES[i]↑(j-1);
for i←1 step 1 until M do
F[i,1] ← OBS[j];
err ← lsolve(A,H,F,N,M,1);
The polynomial might then be evaluated for any new time by the following
subroutine:
Real procedure EVAL(real tim);
begin
real tmp;
tmp←H[N];
for i←N-1 step -1 until 1 do
tmp ← H[i] + tim*tmp;
return(tmp);
end;
Additionally, there are two other subroutines that do a lot of these
things for you: these are LSQ and RSQ. In the previous example, we
could have used LSQ directly with the following call:
err←LSQ(X,OBS,H,M,N);
This takes directly the X array (dimensioned [1:M]) and the OBS array
(also [1:M]) and returns the coefficients in H[0:N]. Again, it returns
zero upon success and a positive number on errors of the above described
flavors.
In the special case that the X[i] happen to be equally spaced, such as
points at equal intervals of time, then you can use another routine, RSQ,
that is even simpler:
err←RSQ(OBS,H,M,N,startingX,difference);
Here we have made the assumption that
X[i] = startingX + i*difference;
So that
X[1] = startingX
X[2] = startingX + difference;
X[3] = startingX + 2*difference;
. . .
X[M] = startingX + (M-1)*difference;
This is very handy for sampled-data functions, where the abcissa are
uniformly placed.
If you decide that some of the points should be weighted differently
from others, you cannot use this set of routines directly. You will
have to use QSOLVE. The idea here is that we add a weighting function,
Wi, to the error form, yielding the following equation:
M N
E = ↔ W[i]*(OBS[i] - ↔ H[j]*G[j](X[i]) )↑2
i=1 j=1
Note here that we must not confuse the weights Wi with the coefficients
Gj, which are sometimes (as above) called "weights" also. The weighting
indicated by Wi represents something about the "goodness" or the validity
of the i'th observation. In this case, when we take the derivative of
the error with respect to the k'th coefficient, Hk, we get the following
set of linear equations to solve:
M N
0 = ↔ W[i]*G[k](X[i])*(OBS[i] - ↔ H[j]*G[j](X[i]) )
i=1 j=1
for k=1, 2, 3, . . . , N
So you can load up the matrices and solve them with the following
SAIL program fragment:
for k←1 step 1 until N do
begin
F[k]←0;
for j←1 step 1 until N do A[k,j]←0;
for i←1 step 1 until M do
begin
F[k]←F[k] + W[i]*OBS[i]*X[i]↑(k-1);
for j←1 step 1 until N do
A[k,j]←A[k,j] + W[i]*X[i]↑(j+k-2);
end;
end;
err←qsolve(N,A,H,F);
This will not be quite as accurate as using LSOLVE, but it allows the
additional freedom of applying weighting functions to data.
Note that there will always be some error in these calculations. Like
for instance, if the data is exactly a function 1.7*X↑3+0.5, a fit with
a fourth-order polynomial may give non-zero coefficients for the other
terms, like the X↑4 term, the X↑2, and the X term. These errors will
depend on the order of the system demanded. In this example, the total
error will not exceed about .01% of any coefficient.
NUMERICAL ANALYSIS - CURVE FITTING
There are various ways of fitting curves to data. Two ways were given
on the previous page, which were Lagrange polynomial and least-squares
polynomials. Here we will examine several more options, probably the
most useful of which is spline curve fitting.
A natural spline of order 3 is a series of piecewise-cubic polynomials
such that at the breakpoints between the polynomials, the values and
the slopes are constant. This is known to be the smoothest curve that
one can fit to the given data that passes through the given
points. It is best used for interpolation of a given curve that only
has known values at certain points. Note that this system is not useful
for extrapolation because the behavior beyond the end points is exactly a
cubic polynomial and will go shooting off to infinity in a cubic manner.
To use the routines, we again start with an array of observations,
OBS[N1:N2] (we usually take N1=0 and N2=N), and times of observations
TIMES[N1:N2]. We usually assume that TIME[i]>TIME[j] if i>j. That is,
that the times are in increasing order. It returns you three arrays that
we call for simplicity B, C, and D, all having the same dimensions as
OBS and TIMES. You can then calculate the interpolated value at any new
time TIM by first finding the value of i for which TIME[i+1]>TIM≥TIME[i]
and computing the following function:
IVAL ← (( D[i]*(TIM-TIMES[i]) + C[i])*(TIM-TIMES[i])
+ B[i])*(TIM-TIMES[i]) + OBS[i];
Be sure to add in the OBS[i] term. The routine that calculates the
values of B, C, and D for you is this:
External procedure CUBNATSPLINE(integer n1, n2;
real array TIMES, OBS, B, C, D);
To help you with the evaluation of interpolated values, the following
routine will find the index i corresponding to the desired value of time:
External integer procedure FIND_SPLINE(real TIM;
real array TIMES; integer n);
So you could make an evaluation routine that looks like the following:
real procedure INTERP(real TIM);
begin
integer i;
real del;
i←find_spline(tim,times,N);
del←TIM-TIMES[i];
return((( D[i]*del + C[i])*del + B[i])*del + OBS[i]);
end;
For that to work, B, C, D, TIMES and OBS must all be dimensioned [0:N].
NUMERICAL ANALYSIS - OPTIMISATION
We have seen already how linear functional optimisation can be done
using the solution of least-squares formulation by matrix inversion.
The only thing remaining then is non-linear optimisation. This is
a rather hairy topic, and should be approached with trepidation.
The routine we have for doing this is the Marquardt modification to the
Gaussian iteration. You must give the function values and the derivatives
and it will return the unknown parameters of the minimum.
The idea here is that (again) given M observations OBS[0:M-1] and a
function G that depends on the observation times TIMES[0:M-1] and on
an unknown parameter vector H[0:N-1], we will try to minimize the
following error function:
M-1
E = ↔ (OBS[i] - G(H,TIMES[i]) )↑2
i=0
An example of this might be trying to fit the decay of the reverberation
in a room in a particular frequency band by a sum of exponentials. One
way to do this might be to define the function G as follows:
N/2-1
G(H,TIM) = ↔ H[2*i] * e↑(H[2*i+1]*TIM)
i=0
Here, the even numbered H parameters (H[0], H[2], etc) are the amplitudes
of the exponentials and the odd numbered ones (H[1], H[3], etc) are the
inverses of the time constants. For this to work right, OBS ought to really
be the backwards-integrated squared decay curve, but let's just pretend
that this is o.k. for now.
To do the minimization at all, we need the derivatives of the function
G also with respect to the parameters H. For the function we have defined
above, these derivatives can be computed as the following:
DERIVS[i,j] = e↑(H[i]*TIMES[j]) for i odd
DERIVS[i,j] = H[i-1] * TIMES[j] * e↑(H[i]*TIMES[j]) for i even
This gives us a matrix of derivatives with respect to each of the variables
we are adjusting, evaluated at each of the times of the observations.
The way this all works is as follows. You supply two subroutines, FUNC
and DERIV. FUNC evaluates the function G at the observation times, and
DERIV evaluates the first derivatives of the function with respect to each
the variables, evaluated at each of the observation times. You then call
the following routine:
External integer procedure MQ(real array H; integer N;
real array G,OBS; integer M;
real array DERIVS);
Where H[0:N-1] is the initial parameter settings (you are strongly advised
to insert some good guesses here), G[0:M-1] is where the program can
store transient values of the function evaluations, OBS[0:M-1] is the
array of observations, and DERIVS[0:N-1,0:M-1] is the matrix of derivatives
of the functions. G and DERIVS do not have to be set at the call, because
the MQ program will call your FUNCT and DERIV to fill them up. It assumes
the following calling sequences for them:
Internal procedure FUNC(integer H,N,G,M);
Internal procedure DERIV(integer H,N;
reference real array DERIVS; integer M);
Note that H and G are given here as integers because MQ is really deep
down inside a FORTRAN program and the array format is not compatible.
This means that, for instance, parameter H[i] must be accessed by the SAIL
expression MEMORY[H+i,real]. This can be used anywhere that you need to
use H[i]. Likewise, you must use MEMORY[G+i,real] to access the array of
function values. So, for our example as given above, the FUNC and DERIV
routines might be the following, assuming TIMES is a global array:
Internal procedure FUNC(integer H,N,G,M);
begin
integer i,j;
for j←0 step 1 until M-1 do
begin
real tmp;
tmp←0;
for i←0 step 2 until N-2 do
tmp ← tmp + memory[H+i,real]*
2.71828183↑(memory[H+i+1,real]*TIMES[j]);
memory[G+j,real]←tmp;
end;
end;
Internal procedure DERIV(integer H,N;
reference real array DERIVS; integer M);
begin
integer i,j;
for j←0 step 1 until M-1 do
for i←0 step 2 until N-2 do
begin
DERIVS[i,j] ← 2.71828183↑(memory[H+i,real]*TIMES[j]);
DERIVS[i+1,j] ←
memory[H+i,real]*TIMES[j]*
2.71828183↑(memory[H+i+1,real]*TIMES[j]);
end;
end;
Luckily, the program MQ uses the original DERIVS array that you supply in
the call to MQ, so that it is a real SAIL array, and can thus be accessed
by the standard array mechanism.
You can then call MQ and it will grind away. One thing you might want to
do is to put print statements into FUNC so that it will tell you what
the error is at each point so you can watch and see if it seems to be
converging. If it diverges and starts getting FLOATING OVERFLOWs, then
you probably have a problem.
NUMERICAL ANALYSIS - MISCELANEOUS
FACTORING POLYNOMIALS
For polynomials A[0:N] of Nth order such that the coefficient A[0]
represents the constant term, A[1] represents the 1st-order term, and
so on, the following routine will factor the polynomial of order up to
99:
External integer procedure POLY(real array A;
reference real array ZR,ZI;
integer N);
It returns as an integer the number of roots actually found. If this
number is negative, it means that you gave the leading coefficient (A[N])
equal to zero (that is a no-no). It returns you that many roots, up to N.
The real parts will be in ZR[1:N] and the imaginary parts will be in
ZI[1:N]. For complex conjugate roots, the imaginary parts will be exact
additive inverses of each other. If the polynomial has multiple roots or
is otherwise ill-conditioned, the program may not be able to obtain all
the roots. If not, the integer returned gives the number of roots that it
did find.
As an example, the following little program fragment
(missing the declarations) reads in from the user the coefficients of
a polynomial and prints out the results:
while true do
begin "LP"
inint(N,"Ordre du polynom");
for i←0 step 1 until N do
inreal(A[i],"Coefficient de X↑"&cvs(i));
i←poly(A,Zr,Zi,N);
outstr(cvs(i)&" racines trouve"&crlf&
" REAL IMAGINAIRE"&crlf);
for j←1 step 1 until i do
outstr(cvg(Zr[j])&" "&cvg(Zi[j])&crlf);
end "LP";
(This service is available as a program called RACINE on the system).
POLYNOMIAL ARITHMETIC
External procedure POLMUL(real array out,in1; integer l1;
real array in2; integer l2);
External procedure POLADD(real array out,in1; integer l1;
real array in2; integer l2);
There are two routines for multiplying and adding polynomials. These
assume that the polynomials are in normal form, that is, a polynomial
of degree N is represented by an array of coefficients, H[0:N], where
the index corresponds to the degree. If you have two polynomials, P1
and P2, of orders L1 and L2, then you can multiply them to produce
a polynomial RESULT of order L1+L2 with the following call:
POLMUL(RESULT,P1,L1,P2,L2);
The RESULT array must be dimensioned [0:L1+L2] (or longer).
Likewise, to add to polynomials to produce a polynomial RESULT of order
L1 max L2, this call is appropriate:
POLADD(RESULT,P1,L1,P2,L2);
Notice that the form is identical to that of POLMUL.
OTHER ROUTINES
There are numerous other ones, such as vector dot products, linear
combinations, Bessel functions, Chi-squared functions, random numbers,
both linear and Gaussian, etc etc etc. Rather than duplicate effort here,
it would probably be best to glance through JAMLIB.SAI[SUB,SYS] and see if
there is anything in there you might need.
SIGNAL PROCESSING
There are a few main branches of signal processing available:
1) Fast Fourier Transforms
2) Filtering
3) Linear Prediction
4) Sound file Input/Output
SIGNAL PROCESSING - FAST FOURIER TRANSFORMS
There are two principal routines for complex FFTs, FRXFM and
FFTARB. The first one is used for data lengths that are strictly powers of
2, and the second one will do any length data, with the provision that it
will run a bit slow for prime lengths. The first routine (FRXFM) is
hand-coded in a most abstruse manner in MACRO, and has an incompatible
calling sequence as follows:
FRXFM(LOG2N,REAL,IMAG);
This must be done in START_CODE in SAIL as follows:
start_code
label lp,rp,ip;
integer save12,save16,save17;
move 1,real; Comment pick up address of 1st word of data;
hrrm 1,rp; Comment from array pointer word;
move 1,imag;
hrrm 1,ip;
movei 1,log2N; Comment get address of length word (integer);
hrrm 1,lp;
movem 12,save12; Comment Save key registers for SAIL;
movem 16,save16;
movem 17,save17;
jsa 16,frxfm;
lp: jump 0; Comment - Address of log2 of number of pts;
rp: jump 0; Comment - Address of 1st wd of real array;
ip: jump 0; Comment - Address of 1st wd of imag array;
move 17,save17; Comment Put the world back;
move 16,save16;
move 12,save12;
end;
After the JSA 16,FRXFM should come the address of the log (base 2) of the
number of words of the transform, then the address of the first data word
of the real array, then the address of the first data word of the imaginary
array. Both these arrays will be overwritten with the complex transform.
The log (base 2) is really simpler than you may think: the log base 2 of
1024 is 10, of 4096 is 12, etc.
Notice that this routine will also take the inverse transform. Before
calling FRXFM, merely set IMAG[i]=-IMAG[i]. That is, just invert the
signs of the imaginary components, then call FRXFM normally.
Luckily for you, there is a simpler way of doing this in the (normal)
case that the data is real. This is done either using RTRAN if the
array is small enough (smaller than 65536 words or thereabouts), or
XRTRAN if you need another factor of two.
External procedure RTRAN(real array X,A,B;
integer log2N; boolean inverse);
If your array is suitably small, RTRAN can be used. If the number of
points is N such that N←2↑log2N, then the arrays should be
dimensioned as follows:
real array X[1:2*N],A,B[0:N];
Comment **** Load X array here ****;
rtran(X,A,B,log2N,false);
On return, the A and B arrays will be filled with good data. Each point
in A or B will be the Fourier Cosine and Sine series. If your sampling
rate, for instance, is F hertz, point A[i] and B[i] will represent
the real and imaginary parts of the transform at a frequency corresponding
to i*F/(2*N). You will immediately notice that this only covers the
frequencies up to half the sampling rate. This is because the input data
are real, meaning that the real part of the transform is even and the
imaginary part of the transform is odd, so that A[i]=A[2*N-i] and
B[i]=-B[2*N-i] (as if those higher parts of the array existed). This
is where the economy of RTRAN comes from: the information over the full
range of 0 to 2N is redundant, and therefore need not be returned to the
user. You may also wonder about the fact that since A and B are
dimensioned 0:N, that means that there are 2N+2 points returned. This may
seem a bit funny since there were only 2N points as input data. This can
be explained by the fact that B[0] and B[N] will always be zero for real
input data.
If the INVERSE parameter is FALSE, the transform of X is taken and the
resulting Fourier cosine and sine coefficients are placed into the
A and B arrays. If INVERSE is TRUE, the inverse transform is done. The
A and B arrays are taken as input data and they overwrite the X array
as result.
More specifically, the transform (INVERSE=FALSE) may be represented
by the following formulae:
2N-1
A[k] = (1/N) ↔ X[j+1] COS(π*j*k/N)
j=0
2N-1
B[k] = (1/N) ↔ X[j+1] SIN(π*j*k/N)
j=0
For k=0, A[0] is just the sum of the input points (and B[0] is zero).
For A[1], it is the cross-correlation of the input waveform and
a one period of a cosine, sampled by 2N points.
If you absolutely have to do a very long transform, such that you cannot
afford to have the X, A, and B arrays existing all at the same time,
then you must use XRTRAN and FRXFM.
External procedure XRTRAN(real array A,B; integer N; boolean evaluate);
The way this works is that you put X[1] into A[0], X[2] into B[0],
X[3] into A[1], X[4] into B[1], and so on. You then call FRXFM on
these data with A in the place of REAL and B in the place of IMAG. After
this, you call XRTRAN with EVALUATE=FALSE and that you will get
back the real transform in A and B.
Likewise, if you wish to do an inverse transform, you call XRTRAN with
EVALUATE=TRUE, then go through and invert the signs of the elements of
the B array (B[i]←-B[i]), then call FRXFM as described above. Your real
data will be found in alternating arrays, so that
X[2*K+1] = A[K]
X[2*K+2] = B[K]
for K=0, 1, . . . , N-1
TRANSFORMS OF ARBITRARY LENGTH
Of the three routines mentioned above, only XRTRAN is capable of dealing
with any length of data. Both FRXFM and RTRAN have the power-of-two length
restriction built in quite strongly. There are two analogous routines for
data of any length. These are FFTARB and ARTRAN.
External procedure FFTARB(real array A,B integer N,Nv,Ks)
This routine actually does multivariate transforms and the reader is
urged to consult the original article for this routine for the details
about doing multivariate transforms:
"Fast Fourier Transform routine with arbitrary factors" CACM Algorithm
339, Richard C. Singleton, CACM 11, #11, November 1968, pages 776:779.
For just doing ordinary 1-dimensional transforms, N=Nv=Ks. For example,
to do our transform as above, we would make the call:
FFTARB(A,B,N,N,N);
Note that this is a SAIL program and expects the arrays to be dimensioned
A[0:N] and B[0:N]. This routine also can be used for inverse transforms
by inverting the sign of the B array.
For doing arbitrary length real transforms, the routine ARTRAN is
available as follows:
External procedure ARTRAN(real array X,A,B; integer N; boolean inverse);
Here, X is dimensioned [1:2*N], and A and B are dimensioned [0:N]. Other
features are identical to RTRAN except that you give the actual number
of points N, rather than its log base 2.
WINDOWING
For whatever reason, people seem to like to apply windows to their
data before transforming it. This can be done with either of two different
programs: WEIGHT and ABWEIG, depending on whether the data is in a linear
array (as in X above) or alternates between two arrays (like in the A and
B arrays before the call on XRTRAN).
External procedure WEIGHT(real array X; integer N;
real array D; integer M);
This computes the "sum of cosines" windowing family that includes the
well-known Hanning and Hamming windows. The weighting function may
be computed as follows:
M
W[i] ← 1 + ↔ D[j] COS(2πij/N) i = 0, 1, 2, . . . N-1
j=1
Some popular windows and their specs:
"Hanning" window: M←1, D[1] ← -1.0
"Hamming" window: M←1, D[1] ← -0.851851851851
2nd order window: M←2, D[1] ← -1.19685
D[2] ← 0.19685
3rd order window: M←3, D[1] ← -1.43596
D[2] ← 0.497537
D[3] ← -0.0615762
4th order window: M←4, D[1] ← -1.566272
D[2] ← 0.725448
D[3] ← -0.180645
D[4] ← 0.0179211
So, for instance, to weight and transform some data of length 2N where
N←2↑log2N, you might do the following:
preload_with -0.851851851851;
real array D[1:1]; Comment Load a Hamming window;
real array X[1:2*N],A,B[0:N];
Comment *** Load up the X array here ***;
weight(X,2*N,D,1); Comment Weight the sequence;
rtran(X,A,B,log2N,false); Comment Transform it;
Notice that autocorrelation and cepstrum may be done by two consecutive
FFTs. Not only that, there is a routine that does it for you:
External procedure AUTCEP(real array X, A;
integer log2N; boolean cepstrum);
Dimensions are X[1:N], A[0:N-1] where N=2↑log2N
Does "windowed" autocorrelation or cepstrum of X and puts it in A. Does
it all by FFT. X and A can be the same array
To be more precise, the autocorrelation (CEPSTRUM=FALSE) is equivalent
to the following sum:
N-k
A[k] = sum X[n]*X[n+k] k=0, 1, 2, . . . N-1
n=1
The cepstrum is the inverse discrete Fourier transform of the log
magnitude spectrum of the points in X. We shall not try to describe
what use it is here, but instead refer the interested user to the
literature (Oppenheim & Schaefer "Digital Signal Processing,"
Prentice-Hall).
SIGNAL PROCESSING - FILTERING
There are three kinds of filtering we usually do with these routines:
1) IIR filtering
2) FIR filtering by convolution
3) FIR filtering by FFT
Since the first one of these is the easiest to explain, I will start
there.
In IIR filtering, there are basically three routines: one for designing
filter coefficients for the standard forms of filters, like lo-pass,
hi-pass, bandpass and bandstop filters, one routine for realizing
these filters, and another routine for doing this same sort of filtering
but with a continuously changing center frequency.
External boolean procedure GFINI(real array H; real Wc1,Wc2,clock;
integer type,class,order);
The filter is specified entirely by the band edge frequencies (Wc1 and Wc2
- latter not used for high and low pass), the sampling rate (CLOCK), the
type of approximation used (Butterworth or Chebychev), and the filter
class (hi-pass, lo-pass, etc), and the order of the filter. This order is
doubled for double-edge filters (bandpass and bandstop). It loads up the
array H with the coefficients with the pass-band amplitude scaled to 1.0.
The array H should be dimensioned to be 1+13*N where N is the order of the
filter. Wc1 and Wc2 are the 3dB points of the filter. Wc2 is not used
for low and high-pass filters. Note that these are in randians per second,
such that Wc=2πf where f is in Hertz.
The types are decoded as follows:
Type=0 Butterworth
Type=1 .5 dB Chebychev
Type=2 1.0 dB Chebychev
Type=3 2.0 dB Chebychev
Type=4 3.0 dB Chebychev
And the classes are decoded thusly:
class=0 Low pass
class=1 High pass
class=2 Band pass
class=3 Band stop
The H array is a mixed mode array and is loaded up in terms of
second-order sections. We will let I stand for the section number
(1, 2, 3, etc). The maximum section number will be half the order
rounded up to the nearest integer, except for bandpass and bandstop,
where it will equal the order.
Each section has 5 coefficients. The output for point N (Y[N]) is
computed from the input points (X[N]) as follows:
Y[N] ← C0*X[N] + C1*X[N-1] + C2*X[N-2] + D1*Y[N-1] + D2*Y[N-2]
The H array is then interpreted as follows:
WORD INTERPRETATION
--------------------------------------------------------------
H[1] (Integer) Total number of sections
H[13*i+0] (Boolean) Order of this section.
0 = FALSE = 1st order, -1 = TRUE = 2nd order
H[13*i+1] (Boolean) Tells if C1=0 or not
0 = FALSE = C1 is nonzero, -1 = TRUE = C1 is zero
H[13*i+2] (Real) C0 - Coefficient of X[N]
H[13*i+3] (Real) C1 - Coefficient of X[N-1]
H[13*i+4] (Real) C2 - Coefficient of X[N-2]
H[13*i+5] (Real) D1 - Coefficient of Y[N-1]
H[13*i+6] (Real) D2 - Coefficient of Y[N-2]
H[13*i+7] to H[13*i+12]
Set to zero by GFINI, used by GFILT for history terms
of X and Y
Once your filter is designed, you can apply it to a sampled-data signal
with GFILT:
External procedure GFILT(real array X,Y,H; integer N; boolean backwards);
This filters the X array (dimensioned [0:N-1]) and places the result
in the Y array (same dimensions). The two arrays can be the same. The
filter coefficients are taken from the H array which is assumed to have
been set up by GFINI. The BACKWARDS flag is normally FALSE, but if TRUE
means to filter the array starting at X[N-1] and going to X[0], i.e.,
in reverse order.
All the history terms are kept in the H array, so that a sampled-data
signal of length longer than N can be filtered in windows of N points
at a time. Any call on GFINI will cause these history terms to be set
to zero.
FIR FILTERING
The first thing needed is how to design FIR filters. To aid this, routines
are provided for low-pass, differentiator, and hilbert transform filters
by windowing. The window used is the Kaiser window, which has a parameter
BETA that determines the passband and stopband ripple. Rather than having
to make this number up out of thin air, there is a routine that will help
you determine both the length of your filter (number of points in the
impulse response) and the value of Beta for these specifications.
External integer procedure DESLPF(reference real beta; real tw,delta);
To use this routine, you must give it the transition width (TW) and the
allowable ripple (DELTA). For instance, if you will tolerate ripple
that is 60 dB down, you give 0.001. Likewise, 40 dB down gives 0.01.
This ripple will be present in both the stop band and in the pass
band. The transition width is in terms of normalized frequency. For
instance, if you are designing a low pass filter that you want to
pass frequencies up to F1, then reject frequencies from F2 on up,
you would give a transition width equal to (F2-F1)/Fs where Fs is
the sampling rate. You cannot give zero here. You must allow some
"slack" space for the filter to change. If you give zero, it will tell
you that the filter you requested will have to be of infinite length
to match up to those specifications. For the differentiator,
the transition occurs near the sampling rate, so that the transition width
is then (Fs/2 - F)/Fs = .5-F/Fs where F is the highest frequency where
accurate adherence to the "ideal" differentiator is required. For the
Hilbert transformer, the transition occurs at zero frequency and at
half the sampling rate, so there are, in fact, two transition regions.
This routine returns you two numbers: Beta and N. Beta is the number
that controls the Kaiser window and 2*N+1 will be the length (in samples)
of the impulse response of the filter. Armed with these two numbers, you
can then call the routines that actually design the filters:
real array Imp[0:N];
External procedure FIRLPF(integer Imp,N; real cutoff,clock,beta);
External procedure FIRDIF(integer Imp,N; real cutoff,clock,beta);
External procedure FIRHIL(integer Imp,N; real cutoff,clock,beta);
These programs return you one "wing" of the impulse response in IMP. The
cutoff frequency is in Hertz, as is the sampling frequency, CLOCK. Beta
and N, you have presumably obtained from DESLPF above. The reason that
only one wing is needed is because for low pass filters, the impulse
response is symmetric, so that IMP[-I]=IMP[I]. Similarly, differentiators
and hilbert transformers are anti-symmetric, so that IMP[-I]=-IMP[I].
There are 4 different routines for this (before you go to an FFT filter)
for these different purposes. These are as follows:
External procedure INTERP(real array X,Y,Imp; integer Np,R,N;
boolean symmetric);
Here, the arrays should be dimensioned X[1-N:Np+N], Y[1:R*Np], and
IMP[0:R*N]. This routine is for doing digital interpolation. That is,
for every point in the input array, you want to produce R points in the
output array. This means you must come up with R-1 new points. The accepted
way of doing this is to use a low-pass filter for this interpolation.
You can just call FIRLPF above to design one of the appropriate length
and use INTERP to do the interpolation for you. We have rounded the
filter length up to the nearest multiple of R because it doesn't take
any more computing to do so. For interpolation, the boolean SYMMETRIC
should be set to TRUE. This refers to the impulse response. Note that
the program can be run with R=1, which would be just simple FIR filtering
with no interpolation. In this case, you might want to do something
besides symmetric filters. If you set SYMMETRIC to FALSE, it takes the
filter to be anti-symmetric, as is needed for a band-limited differentiator
or Hilbert transformer. Do note that the length of the wing of the impulse
response IMP is of length R*N, not simply of length N.
External procedure RESAMP(real array X,Y,Imp; integer Np,R,N);
Dimensions are X[1-R*N:R*Np+R*N], Y[1:Np], and IMP[0:R*N]. This routine
is normally for time decimation (i.e., reducing the sampling rate). It
is sort of the opposite of INTERP above. That is to say, that for every
R points in the input array, it produces one point in the output array.
As before, the length of the impulse response in IMP has been rounded
up to the nearest multiple of R. The idea here is that to do decimation
in time (sometimes called "resampling" or "downsampling"), you generally
want to filter the signal with a low-pass filter to eliminate all frequencies
above .5/R times the original sampling rate to prevent aliasing. To do
this, you can use FIRLPF to design the impulse response of a filter with
a cutoff of somewhat less than .5/R times the sampling rate. In this routine,
R must be strictly greater than unity, and the filter is forced to be
symmetric.
External procedure HILBER(real array X,Y,Imp; integer Np,N);
Dimensions are X[1-N:Np+N],Y[1:Np],Imp[0:N]. This assumes that the
filter in IMP is a Hilbert transformer. In this case, it will be
anti-symmetric and the even numbered points will all be zero, allowing
the routine to economize the number of multiplications.
External Integer procedure SRFILT(real array In; integer sip;
real array out; integer sop,Nop,num,den;
real array Imp; integer Np);
This routine is sort of a combination of INTERP and RESAMP. It is
normally for sampling rate conversion by a factor of NUM/DEN.
**** Look at SRCONV to figure this thing out ****
Routine to help do sample rate conversion by a factor of NUM/DEN.
Takes input array IN and filters it by the impulse response IMP[0:Np-1]
(assumed to be one side of a symmetric, odd-length impulse response)
and puts the result in OUT. The sample number of the first element of
the input array is SIP, of the output array is SOP. There are NOP
output points. This program inner loop uses SRFILT to do this:
⊃ Nip←Nop*num/den+21; ⊃ constant depends on filter;
⊃ for sop←0 step Nop until oplen do
⊃ begin "RI"
⊃ sip←sop*den/num-10; ⊃ That constant depends on the filter used;
⊃ sndin(chan,sip,Nip,IH,X,iflag);
⊃ srfilt(X,sip,Y,sop,nop,num,den,imp,Np);
⊃ sndout(ochan,Nop,OH,Y,oflag); ⊃ In/out arrays must be different;
⊃ end "RI";
⊃
⊃ This routine would be equivalent to the following SAIL program:
⊃
⊃ for i←0 step 1 until Nop-1 do
⊃ begin
⊃ off←((i+sop)*den) mod num;
⊃ stp←((i+sop)*den) div num;
⊃ npts←(Np-off-1) div num;
⊃ for j←0 step 1 until npts do
⊃ Out[i]←Out[i]+Imp[off+j*num]*In[stp-j-sip];
⊃ off←-off;
⊃ npts←(Np-off-1) div num;
⊃ for j←1 step 1 until npts do
⊃ Out[i]←Out[i]+Imp[off+j*num]*In[stp+j-sip];
⊃ end;
All these filters will be exactly linear phase.
FIR FILTERING BY FFT
External procedure FFTINI(Real array Imp; integer Nfp; integer symmetric);
External procedure FFTFLT(Real array X,Y; integer Np);
External procedure FFTREL;
Filtering by FFT is much the same as with any other FIR filtering program,
except that it is broken up into three parts: the first routine, FFTINI,
takes the FFT of your impulse response and stores it internally in arrays
that it allocates itself. The next routine, FFTFLT, actually does the
filtering and should be applied as many times as necessary to exhaust the
input data. For each Np points in the input array, Np points are produced
in the output array. To account for the transient response of the filter,
then, you must skip the first and last Nfp points of the output array
(except in the case of a one-sided impulse response, where you skip only
the first Nfp points of the output array: see below).
When the input data is exhausted, FFTREL should be called to release the
internal storage that was used to hold the transform of the impulse
response. Notice that this means that you can only run one filter at a
time in that there is only a unique transformed impulse response in the
program. If you call FFTINI a second time, it will clobber these data.
SYMMETRIC is a 3-valued number here: +1 for symmetric data, such that
IMP[j]=IMP[-j]. -1 is for antisymmetric impulse responses, such that
IMP[j]=-IMP[-j]. SYMMETRIC=0 means that the impulse response has no
symmetry and the entire response is contained in IMP[0:Nfp-1]. This last
case is for doing things like convolving a sound source with the impulse
response of a concert hall. Note that if the impulse response has any
symmetry (symmetric or anti-symmetric), there will be Nfp points at
the beginning of the output array and Nfp points at the end of the output
array that will not be valid, because it will consider the points outside
the range of the input array to have been zero. Likewise, if the impulse
response is one-sided (has no symmetry, SYMMETRIC=0), then there will
be Nfp points at the beginning of the output array that will be as if
the Nfp points before the beginning of the input array were used in the
computation but were zero.
SIGNAL PROCESSING - LINEAR PREDICTION
The linear prediction routines can be divided into roughly three different
groups:
1) Routines to compute the filter coefficients.
2) Routines to apply the filter or inverse filter.
3) Routines to manipulate the coefficients.
FILTER COMPUTATION
There are three kinds of linear prediction supported in this package:
(1) The Wiener-Levinson recursion ("autocorrelation" method) (2) The
covariance method (Yule-Walker equations) (3) The Burg method ("Maximum
Entropy" method). We will take these one at a time:
External procedure AUTC(real array X; integer N; real array R; integer M);
External real procedure AUTSOL(real array R,A,K; integer M);
The process of linear prediction filter computation by the autocorrelation
method is to first compute the autocorrelation coefficients of the (windowed)
sound file, then to solve the Toeplitz form simultaneous linear equations for
the filter coefficients. The program for computing the autocorrelation
coefficients is AUTC, and it computes the following function:
N-i
R[i] = ↔ X[j]*X[j+i] for i=0, 1, 2, . . . , M
j=1
The dimensions are R[0:M] and X[1:N]. For normal use, the signal, X,
should have been windowed using, for example, a Hamming window. The
subroutine WEIGHT described previously can do this.
The routine AUTSOL does the solution to the normal equations by the
Levinson recursion. It returns the filter coefficients in either
direct form (A[0:M]) or as reflection coefficients (K[1:M]). The
direct form is useful is you wish to factor the polynomial to observe
the roots, and the reflection coefficients are useful because they
result in the filtering operation with the minimum roundoff error.
Note that A[0] will always be set to 1.0 by the solution routines.
Note also that the Levinson recursion is not numerically stable and
can give bad results for higher orders, or if the eigenvalue spread
of the Toeplitz matrix is too high.
Note that the AUTSOL routine can be used for two other purposes:
(1) to modify a given set of K-parameters that have already been
determined previously (2) to design a filter (that is, a set of
K-parameters) that realizes a certain spectrum. To modify a set
of K-parameters, you first convert them to direct form (or A-parameters)
using the routine K2A (described below), then transform them to the
spectral domain using any of the FFT routines described earlier. The
resulting spectrum is then the inverse of the spectrum of the filter
in recursive form. You may then modify this spectrum as much as you
want, then calculate the autocorrelation function from this spectrum
by taking the magnitude squared, then the inverse FFT. This set of
autocorrelation coefficients may then be fed directly to AUTSOL to
compute a new set of K-parameters for the modified spectrum. Likewise,
we do not need to start with a spectrum from an existing set of
K-parameters, we can start immediately from a spectrum (magnitude
squared) that you draw by hand, or invent in any other way. The
inverse FFT of this spectrum thus gives autocorrelation coefficients
that can be used via AUTSOL to produce the filter coefficients.
The program KWARP does exactly this process to modify or create filters.
To use the covariance method, we simply calculate the covariance matrix,
then set up the normal equations (sometimes called the Yule-Walker
equations). The solution is then obtained by matrix inversion, as was
described before. The procedure to calculate the covariance matrix is
as follows:
External procedure COVAR(real array PHI,DUMMY; integer M,N,X);
PHI should be dimensioned [0:M,0:M], DUMMY[1:2*M+2], X[0:N-1]. This
computes the following sum:
N-1
PHI[i,j] = ↔ X[n-i]*X[n-j] for 0 ≤ i,j ≤ M
n=0
The normal equations, however, do not use the 0th row and column,
so we must reformat this into an M*M matrix and a 1:M column matrix.
This can be done as follows:
real array P[1:M,1:M],B[1:M];
for i←1 step 1 until M do
begin
B[i]←PHI[0,i];
for j←1 step 1 until M do
P[i,j]←PHI[i,j];
end;
The matrix may be inverted by using QSOLVE as follows:
error ← QSOLVE(M,P,Ap,B);
This gives you the predictor coefficients, Ap[1:M], but not the reflection
coefficients. Note that Ap[0] is not computed, but is universally 1.0.
The filter is not necessarily stable.
The third method of computing the predictor coefficients is by Makhoul's
modification of Burg's method. This is like the covariance method, except
that the stability of the filter is guaranteed. The calling sequence
is similar, in that you use COVAR as described above, but the solution
for the filter is done by another routine:
External procedure CLATT(real array PHI,A,KP; Integer P);
Dimensions are PHI[0:P,0:P], A[0:P], Kp[1:P]. This solves directly
for the filter coefficients in direct form (A[i], with A[0]=1.0),
and for the reflection coefficients (Kp[i]).
FILTER APPLICATION
There are various forms of filtering routines, depending on exactly what
you want to do with it.
External procedure LATICE(real array X,Y; integer N;
real array K,H; integer M; boolean inverse);
This is the straightforward program for filtering, by the K-parameters,
a segment of a signal. It filters the data in X[1:N] by the filter whose
reflections coefficients are in K[1:M], and places the result in Y[1:N].
(X and Y may be the same array]. The history array, H[1:M], should be
set to zero before the first call. If the data is longer than will fit
into a single array, it may be divided into units of N points and the
history array, H, should not be disturbed between successive calls. It
preserves the state of the filter after the pass through the array.
If INVERSE is TRUE, the inverse filter will be used, thus (probably)
whitening the signal in X. If INVERSE is FALSE, the "forward" filter
is used, and X should probably have a fairly flat spectrum to start with.
In the case that the filter is not constant, but changes with time, there
is another procedure for this case:
External procedure VLATT(real array X,Y; integer N;
real array K; integer koffset; real array H;
integer array cmprsses; integer M,Nks; boolean inverse);
Filters from the X array to the Y array, both of which must be N
points long. Uses the reflection coefficients from K[1:M,1:NKS].
The columns (K[*,i]) represent one filter each, of length M. This
is like the K[1:M] used above in LATICE, except that there are
more than one of them. Each filter is assumed to extend over
CMPRSSES[i] points. We will start KOFFSET points into the first
filter. KOFFSET=0 means you start exactly at the first filter.
The K parameters are continuously interpolated between adjacent filters.
All other parameters are as before. The program stops when the N input
points are exhausted, and NKS better be high enough so that there are
sufficient filters to do the entire range, in that no check is done
inside the subroutine to see if we have run off the end of either the
K or the CMPRSSES array.
FILTER MANIPULATION
The things one might attempt to do with are filter (that we will discuss
here) are as follows:
1) Factoring a filter to find the roots
2) Getting the spectrum of a filter
3) Placing it into parallel or cascade form
For all of these, the filter must be in the direct form (denoted by
the A[0:M] array above). If you only have the K parameters, you must
convert them to direct form by use of the following routine:
External real procedure K2A(real array Kp,Ap;
integer M; boolean inverse);
The dimensions are Kp[1:M] and Ap[0:M]. If INVERSE is FALSE, the routine
converts the reflection coefficients to the direct form coefficients.
(Note that Ap[0] will always be 1.0). If INVERSE is TRUE, it converts
the direct form coefficients to the reflection coefficients, but it is
a numerically unstable transform, so that only very low orders can be
inverse transformed. The forward transform (INVERSE=FALSE) is numerically
perfectly stable, even up to high orders.
You can now, for instance, use POLY (described earlier) to factor the
polynomial into its component roots. The only thing to remember is that
the filter represented by the A[0:M] coefficients is in terms of negative
powers of Z. i.e., to filter something by A, you would have to do the
following:
M
Y[i] = ↔ A[j]*X[i-j] INVERSE filter
j=0
M
Y[i] = X[i] - ↔ A[j]*Y[i-j] ALL-POLE filter
j=1
This means that if you just factor the polynomial, you will get the
inverses of the roots. To get the correct sense, you can either invert
the roots after factoring, or reverse the order of the coefficients.
The following program fragment shows the conversion of a row of K coefficients
to A coefficients, then the reversal of the ordering and subsequent
factorization:
k2a(Kr,Ap,M,false);
for i←0 step 1 until M do
tA[i]←Ap[M-i];
j←poly(tA,Zr,Zi,M);
Likewise, to take the spectrum of the filter, we might do the
following (assuming "!" is a comment character):
for log2N←6 step 1 until 13 do
if (1 ash log2N)≥M then done; ! Find nearest power of 2;
N←1 ash log2N; ! Number of points represented;
begin
real array Ax[0:2*N-1],Rl,Im[0:N];
k2a(Kp,Ap,M,false); ! May not be needed;
arrclr(Ax); ! Is redundant in this case;
arrblt(Ax[0],Ap[0],M+1); ! Not necc. if Ap array is big;
rtran(Ax,Rl,Im,log2N,false); ! Take FFT;
! At this point, we have INVERSE spectrum in Rl-Im ;
for j←0 step 1 until N do
Im[j]←1/(Rl[j]↑2+Im[j]↑2);
! Change to magnitude, take inverse;
Many of these steps may not be necessary in your particular case. You
must remember, however, that just transforming the coefficients will
yield the INVERSE filter, and that for the all-pole model, you must
invert this spectrum.
To convert the filter to cascade form, you must first factor the polynomial
as given above, then reassemble the coefficients into second-order
sections by locating the complex conjugate roots. This may be done
by the following subroutine:
Procedure COMBINE(Real array Rls,Ims,Ai,Bi; integer M);
begin
boolean array taken[1:M];
integer ind,upto,jj;
arrclr(taken);
ind←1;
for upto←1 step 1 until M do
begin "FFCR"
if taken[upto] then continue "FFCR";
if abs(Ims[upto])<.0000001 then
begin "ISR"
for jj←upto+1 step 1 until M do
if abs(Ims[jj])<.0000001 then done;
if jj≤M then
begin "ISANO"
Ai[ind]←-Rls[upto]-Rls[jj];
Bi[ind]←Rls[upto]*Rls[jj];
taken[jj]←true;
end "ISANO"
else begin "NOTAN"
Ai[ind]←-Rls[upto];
Bi[ind]←0;
end "NOTAN";
end "ISR"
else begin "ISI"
for jj←upto+1 step 1 until M do
if abs(Ims[upto]+Ims[jj])<.0000001
∧ abs(Rls[upto]-Rls[jj])<.0000001 then done;
if jj>M then usererr(0,0,"Complex filter?????");
Ai[ind]←-2*Rls[upto];
Bi[ind]←Rls[upto]↑2+Ims[upto]↑2;
taken[jj]←true;
end "ISI";
ind←ind+1;
end "FFCR";
end;
This routine takes the real and imaginary parts of the roots in
Rl[1:M] and Im[1:M] and it produces the parameters for second
order sections, where each section is then of this form:
1
Hi(z) = ---------------------
-1 -2
1 + Ai Z + Bi Z
(This is the all-pole form that is shown. The inverse filter is just
exactly the inverse of this).
The total filter then is the product of these sections, H1*H2* . . .
Additionally, one can then take this form, and transform it further
to parallel form by use of the Heavyside expansion. The following
routine performs this, assuming that there are no duplicated roots
(i.e., the multiplicity of each root is exactly one):
External integer procedure PFRACT(real array ai,bi; integer M;
reference real array ci,di);
This returns the parameters Ci and Di such that the total filter
is now the following:
-1
M Ci + Di Z
H(z) = ↔ --------------------
i=1 -1 -2
1 + Ai Z + Bi Z
This is then the parallel form of the filter. If the filter was computed
by linear prediction in any of the forms described above, it is extremely
unlikely that multiple roots will obtain for real-world input signals.
You can force a multiple root by synthesizing the correct input signal,
but in normal circumstances, any noise will cause the roots to split,
thus making this partial fraction expansion possible.
SIGNAL PROCESSING - SOUND FILE IO
Routines: JAMA:SNDIO.FAI[LIB,JAM]
JAMA:RHEAD.SAI[LIB,JAM]
JAMA:WHEAD.FAI[LIB,JAM]
There are four parts to sound file IO: reading sound file headers,
reading the samples, writing sound file headers, writing the
samples. There is a separate subroutine for each of these things.
First off, about sound file headers:
HEADERS
The format of a sound file header is the following:
WORD MEANING VALUES COMMENTS
------------------------------------------------------------------------
0 Magic number 525252525252 This means the file has a header
1 Sampling Rate e.g., 25600 Integer - in Hertz
2 Packing 0, 1, or 3 0 is 12 bit, 1 is 18(16) bit,
3 is 36-bit floating-point
3 Number chans 1, 2, or 4 Mono, Stereo, or Quad
4 Max. amp. ≥0 Floating point number
5:63 Reserved
64:127 Comment Asciz text comment
The 0th word is a magic number that could not be a sample under virtually
any circumstances, so it flags the fact that this sound file has a header.
(N.B. - we may change this in the future to be 525252 in the left half and
an integer in the right half that indicates header version number).
The 1st word is an integer that represents the sampling rate per channel
of the sound file in samples per second per channel. For example, two
common sampling rates are 25600 samples per second and 32000 samples
per second. Again, this is per channel sampling rate, so the actual
total number of samples per second will be this number times the number
of channels. The 2nd word is the packing mode, again an integer. There
are currently three packing modes recognized - 0 for 12-bit samples
packed 3 to a word, 1 for 18-bit samples packed 2 per word, and 3 for
36-bit floating point samples packed 1 to a word. If one has 16-bit
samples, they are packed as the high-order 16-bits of each 18-bit
halfword. The bytes (12, 18, or 36-bit) are packed such that the first
byte is in the high order bits of a word, just like is required for
the PDP-10 byte instructions. The 3rd word of the header is a number of
channels in the sound file. The currently recognized values of this
are 1 (Mono), 2 (Stereo), and 4 (Quadraphonic). It is possible that 8
might by added in the future. In the bulk of the sound file itself,
the samples will be interleaved, that is, first will come sample 1 of
channel 1, then sample 1 of channel 2, then sample 1 of channel 3, and
so on.
The 4th word of the header is a positive floating-point number that represents
the maximum amplitude of the sound file. This is zero if the maximum amplitude
has not been calculated. If the sound file contains floating-point samples,
then this number will be just the maximum of the absolute value of all the
samples. If the sound file is integer (12 or 18 bit 2's complement integers),
they are considered to be fractional numbers, that is, they are considered
to be between -1 and +1-ε. Note that this is at odds with other representations,
such as are found in MUSIC V and MUSIC, where the sample as an integer
number is represented as a fixed-point whole number. For instance, in
MUSIC, 12-bit files have ranges from -2048 to 2047. As far as headers
are concerned, however, and as far as almost all sound-processing programs
are concerned, samples are considered to be between -1 and +1-ε regardless
of the packing. The only exception is 36-bit floating point, which is
capable of holding any number representable on the PDP-10.
We will give here two examples of sound file formats to demonstrate
what they look like in more detail:
Mono file, 18-b samples:
WORD DATUM COMMENT
-----------------------------------------------------------------
0 525252,,525252
1 62000 (Octal for 25600 samples/second)
2 1 18-bit integer packing
3 1 Mono
4 200421,,320713 (Octal for 0.534)
5:63 undefined
64 517215,,771350 Asciz /Short comment/
65 203075,,766732
66 627356,,400000
67:127 0
128 0,,3 Sample 1,,Sample 2
129 5,,8 Sample 3,,Sample 4
130 4,,1 Sample 5,,Sample 6
131 -1,,-4 Sample 7,,Sample 8
. . .
Stereo file, 12-b samples:
WORD DATUM COMMENT
-----------------------------------------------------------------
0 525252,,525252
1 76400 (Octal for 32000 samples/second)
2 0 12-bit integer packing
3 2 Mono
4 200771,,260101 (Octal for 0.987)
5:63 undefined
64 517514,,571312 Asciz /Stereo file/
65 675014,,664730
66 624000,,000000
67:127 0
128 777700057776 Sample 1/chan 1
Sample 1/chan 2
Sample 2/chan 1
129 001177720023 Sample 2/chan2
Sample 3/chan 1
Sample 3/chan 2
. . .
And so on.
Now that you know what a sound file looks like, the next part is how can
you access it? There are four programs for this: SNDIN and SNDOUT for
reading and writing samples, and RHEAD and WHEAD for reading and
writing headers. Let us start with input:
External procedure RHEAD(integer chan;
reference integer iflag,tlen,trap,clock,inchns;
reference real maxamp);
First, you open the disk file in the usual way in mode '17 (binary
dump mode). Then you call RHEAD with the channel number that you
opened it on. If the file has a header, it will print out the header
data and will return these data to you in the variables given. CLOCK
will contain the sampling rate, INCHNS will contain the number of
channels, TLEN will be the length of the file in samples, and
MAXAMP will be the maximum of the absolute values of the samples on a
scale of -1 to +1-ε (except for floating-point files, which may be
any range). TRAP is an error flag and is normally set to FALSE to indicate
that everything is ok. If the header is garbaged, then TRAP will be
set to TRUE and you probably should not try to read this file until
it is fixed. If the file does not have a header, RHEAD will query the
user as to the various parameters (clock rate, number of channels, etc).
The number of samples (TLEN) will be the correct number of samples
regardless of whether the file has a header or not: that is, if there
is a header, the number of samples is the number of words in the file,
minus 128, multiplied by the number of samples per word. If it does not
have a header, the subtraction by 128 is ommitted.
The last unexplained word is IFLAG. This is a micro-coded word for use
with SNDIN and SNDOUT to control the packing and unpacking of the file. It
is decoded as follows:
BIT OCTAL COMMENT
----------------------------------------------------------------------
0 400000,,0 SNDIN returns ON if error, OFF upon success
1 200000,,0 This file has a header
2 100000,,0 This is last output (SNDOUT only)
3:5 070000,,0 Packing mode of the file
0 for 12 bit
1 for 16 bit (halfwords, left adjusted)
2 is reserved
3 for 36-bit floating point
6:8 007000,,0 Data type of X array (Normally = 3)
0 for 12 bit (right adjusted, sign-extended integer (12-bits))
1 for 16 bit (right adjusted, sign-extended halfword (18-bit integer))
2 is reserved
3 for 36-bit floating point
9:10 000600,,0 Number of channels (not currently used, but may be someday)
1 for MONO
2 for STEREO
4 for QUAD
11:35 000177,,777777 Length of file in samples (if known)
RHEAD will set this word up for you all by itself just in a manner that
SNDIN will be happy with. If you want to know the packing of the file,
you can get it with the SAIL call as follows:
pack←ldb(point(3,iflag,5));
which will extract the packing code from the flag word. Normally, you
are set to convert the sound file into 36-bit floating-point for computation,
but if for some reason, you with to keep the samples in unconverted integers,
you can set the byte in bits 6:8 to whatever you want.
To actually read the file, you must call SNDIN with the following
specifics:
External procedure SNDIN(integer channel,sample_number,N;
real array H,X; reference integer iflag);
real array H[1:13];
Where CHANNEL is the IO channel number the file is open on, SAMPLE_NUMBER
is the number of the first sample you wish to read (starting at zero),
N is the total number of samples you wish to read, X is a real array
where the data read in will be placed, and IFLAG is the flag word
described above. If the flag word specifies that the samples are not
to be converted to floating-point, but are instead to remain in
integer, then you should declare X as an integer array. H is a history
array that contains various per-file data, like the current word
we are on, what block number, etc etc. The first word of this should
be set to zero before reading the file (like when RHEAD is called).
For instance, here is a complete procedure (minus declarations) for
opening and reading a sound file:
while true do
begin "SFIL"
outstr("Input SND file: ");
sndfil←default(inchwl,snddev,"SND","DSK");
outstr("Reading from "&snddev&":"&sndfil&crlf);
open(1,snddev,'17,0,0,200,brk,eof);
lookup(1,sndfil,fail);
if ¬fail then done "SFIL";
outstr("File "&snddev&":"&sndfil&" not found."&crlf);
end "SFIL";
rhead(1,iflag,tlen,trap,clock,nchns,maxamp);
iH[1]←0;
Np←2000;
begin "ALRR"
real array X[1:Np];
for bgss←0 step Np until tlen-1 do
begin "RR"
if bgss+Np≥tlen then rNp←tlen-bgss
else rNp←Np;
if rNp<1 then done "RR";
sndin(1,bgss,rNp,iH,X,iflag);
Comment ******* YOUR VERY OWN PROCEDURE GOES HERE. DATA IN X ******;
end "RR";
end "ALRR";
Close(1);
Release(1);
Inside the inner loop at the comment, you may place procedures to
do any kind of processing desired on the sound.
OUTPUT
For sound output, there are the routines WHEAD and SNDOUT. WHEAD does
not really write a header, it merely formats it. You must write it
with an ARRYOUT operation. SNDOUT has the same calling sequence as SNDIN
with a few exceptions: (1) the starting sample number of the transfer
is not included. All outputs are strictly contiguous. Each call on
SNDOUT appends that many more words to the output file (2) the OH history
array is the honest-to-goodness output buffer, so for efficiency, it
should be fairly long (like 2500 words or so). SNDOUT will use however
much of it there is. If you do not give a buffer sufficiently large (at
least 140 words), you will shortly get a monitor error message
ADDRESS CHECK FOR DEVICE DPAx
This error also occurs if the disk is not open in mode '17.
The declaration for the routine is as follows:
External procedure SNDOUT(integer channel,N;
array H,X; reference integer oflag);
Since you don't have RHEAD to set up the flag word for you, you must
set it up yourself. Again, bits 6:8 refer to the data format of the
X array, normally considered to be floating-point, and bits 3:5 refer
to the format of the file that is being written. The calling sequence
for WHEAD is a bit tricky also:
External procedure WHEAD(reference real array head;
integer clock,pack,nchans; real maxamp; string id);
This formats the 128-word buffer in HEAD to be an image of a sound
file header with the given data: CLOCK is sampling rate in samples
per second (per channel), PACK is the sound file packing mode,
NCHANS is the number of channels, MAXAMP is the maximum absolute
value, and ID is a string comment (up to 319 characters in length).
You must write the header onto the file using ARRYOUT before any
sound output is done. If you wish to put the maximum amplitude of
the file into the header also, you must rewrite the header after all
samples have been sent.
To flush out the remaining words in the output buffer, OH, for SNDOUT,
you must make one final call to SNDOUT with bit 2 of the output flag
word turned on. (Don't forget to turn it off before you attempt to
write another file). The following is a template for writing a sound
file:
real array head[1:128];
real array OH[1:2500];
while true do
begin "S2FIL"
outstr("Output SND file: ");
sndfi2←default(inchwl,sndde2,"SND","DSK");
outstr("Writing on "&sndde2&":"&sndfi2&crlf);
open(2,sndde2,'17,0,0,200,brk,eof);
enter(2,sndfi2,fail);
if ¬fail then done "S2FIL";
outstr("Can't ENTER "&sndde2&":"&sndfi2&crlf);
end "S2FIL";
oflag←'013000000000; Comment 18-bit integer file;
pack←1;
maxamp←0;
whead(head,clock,pack,1,maxamp,"Example output - NOT FINISHED YET"&crlf);
oh[1]←0;
arryout(2,head[1],128);
Np←2000;
for bgss←0 step Np until tlen-1 do
begin "ALAR"
real array X[1:Np];
if bgss+Np≥tlen then rNp←tlen-bgss
else rNp←Np;
Comment ****** GENERATE X[1] . . . X[Np] HERE *******;
for i←1 step 1 until rNp do
maxamp←maxamp max abs(X[i]);
sndout(2,rNp,OH,X,oflag);
end "ALAR";
oflag←oflag lor '100000000000;
sndout(2,0,oh,head,oflag); Comment Flush remaining buffer;
whead(head,clock,pack,1,maxamp,"Example output");
useto(2,1); Comment Update header for maximum amplitude;
arryout(2,head[1],128);
close(2);
release(2);
One must be careful that the information given to SNDOUT for the
formatting of the output sound file is consistant with that given to WHEAD
to format the header. It is a quite common error to tell SNDOUT one
packing mode and WHEAD another one. This results in a particular form of
astonishment.
DISPLAY
Routines: JAMA:DPY.FAI[LIB,JAM]
JAMA:FAISUB.FAI[LIB,JAM]
JAMA:PNTPLT.SAI[LIB,JAM]
The display routines at IRCAM are in a somewhat funny state, in that
they work only in conjunction with a PDP-11 which has graphics
capabilities. Additionally, we cannot even access all the features
of the VT-11 graphics processors, in that no provision for light
pen interrupts or any information returning from the PDP-11 is
currently installed. With these caveats, let us proceed:
The graphics routines model the screen as a square with coordinates
in both the X and Y directions from -512 to 511. Nothing reasonable
is guaranteed to happen if you go outside of this range.
The routines work in terms of "pictures." You allocate space for
a picture, then call the display building routines as you want. You
can then post a picture on the screen, or write it out on the
disk (with or without first having posted it on the screen). The picture
can then be erased, or de-allocated. We will describe each of the
routines in order.
To create a picture, you use the following routine:
External procedure DSETUP(integer nwds; reference integer id);
This claims NWDS out of free storage (i.e. - it enlarges your program
by that many computer words) and puts a pointer to this storage
area in the integer variable ID. The variable is called the "picture
identifier," or more briefly, the "picture id." All subsequent reference
to the picture is made through this identifier. You can estimate the
number of words required as the number of points or vectors you wish to
draw. It never hurts to use more than is necessary. The PDP-11 can only
hold about 4000 vectors in its memory. The call to DSETUP must be made
before any other display routines are called. The routine also initializes
the picture to a blank, or empty picture.
To build the picture, there are the following simple routines available:
External procedure RVECT(integer id,dX,dY);
External procedure RIVECT(integer id,dX,dY);
External procedure AVECT(integer id,X,Y);
External procedure AIVECT(integer id,X,Y);
External procedure RPT(integer id,dX,dY);
External procedure APT(integer id,X,Y);
External procedure DTEXT(integer id; string text; real scale,angle);
The simplest way to make a picture is with vectors. There are what
we call "visible" and "invisible" vectors. If we think of the process
of drawing the picture as being analagous to drawing a picture by hand
with a pen, an invisible vector means "raise the pen and move it to
position (X,Y)", and a visible vector means "lower the pen and move it
to position (X,Y)." When a picture is initialized, the "pen position" is
set to (0,0), which is the center of the screen. You may reposition
the "pen" with the AIVECT call as follows:
AIVECT(id,Xc,Yc);
This will move the pen from its current position to (Xc,Yc) without
displaying anything. Likewise, to make a mark on the screen, AVECT
can be called (with the same argument list) to move the pen from
its current position to the new position (Xc,Yc) and leave a straight
line in between.
These coordinates so far are absolute. We also have what are called
"relative" vectors, which move from the current position by a
displacement amount in both X and Y. For example,
RVECT(id,10,5)
will cause the pen to be moved from its current position 10 points to
the right and 5 points up, leaving a straight (slanted) line as it
goes. RVECT(id,-13,45) would move 13 points to the left and 45 points up.
We can draw a square of size 12 at the current position by the following
sequence:
RVECT(id,12,0);
RVECT(id,0,12);
RVECT(id,-12,0);
RVECT(id,0,-12);
This will draw a square, the lower left-hand corner of which will be at
the current position.
RPT and APT are also available to draw single points. They reposition
the pen (either relatively or absolutely) and at the new position, they
draw exactly one dot. This may or may not be very visible on the screen.
All the time you are making these calls, nothing will be shown on the
screen. You have to explicitly post the display on the screen with
the following command:
External procedure WRITE(integer id,pog);
The POG argument is currently inoperative and can be set to zero. At
some point in the future, it may become an "overlay" number, so that you
will be able to update one of several currently displayed pictures, but
that argument is currently ignored.
To draw text (characters) on the picture, the following call is useful:
External procedure DTEXT(integer id; string text; real scale,angle);
It starts at the current pen position. The lower left-hand corner of
the first character will be placed at the current pen position. If SCALE
and ANGLE are both set to zero, the characters used are about 15 by 10
points on the screen (it uses the hardware character generator on the
VT11 display). Angle is denivelation in degrees. For instance, ANGLE
set to 90 causes the text to be vertical, with the baseline on the
right. 180 causes the text to be upside-down, right-to-left. SCALE set
to zero or 15 give "normal" size characters. SCALE set to, for instance,
30 will give double size characters. At all angles and all scales, the
lower left-hand corner of the first character (from its frame of reference)
will be placed at the pen position. After this call, the pen position is
set to draw the next character, if there were one. It does not restore the
pen position to the original point on the screen.
When you picture is complete, you may wish to write it out on the disk.
To do this, you must first open a file for writing in mode '17, then
you call the following routine:
External boolean procedure DWRITE(integer id,chan);
You must CLOSE and RELEASE the channel after this or else the file will
not exist. Once the file exists on the disk, you can redisplay it
with the "D" command in the S program, or you can plot it on the
versatek with a combination of programs, starting with PLOT.
Additionally, after all operations are finished with the picture, you
can release its storage with this routine:
External boolean procedure DRELS(reference integer id);
If you just want to re-initialize the picture without deleting the
storage for it, you can use this routine:
External procedure BUFCLR(integer id,nwds);
The NWDS parameter is the same number that you originally gave to DSETUP.
Since this number is not stored in the picture itself, you must re-tell
the routines how large the picture storage is.
These are all the routines needed to draw pictures on the screen. The
following program fragment draws a picture, posts it on the display,
then writes it out on a disk file. (Assume the "!" is a comment character).
define π="3.14159265";
dsetup(id,500); ! Claim a buffer for the picture;
aivect(id,500,0);
for i←1 step 1 until 300 do
avect(id,500*cos(2*π*i/300),500*sin(2*π*i/300));
write(id,0); ! Post the display on the screen;
open(1,"DSK",'17,0,0,200,brk,eof);
enter(1,"FOO.PLT",fail); ! Default extension is .PLT ;
dwrite(id,1); ! Write picture out;
close(1);
release(1);
drels(id); ! Release storage entirely;
This is enough to get you started. Additionally, you may want to use
some of the more advanced plotting routines to take some of the burdon
of axis generation and such off the poor programmer. For these, we use
two arrays which specify various data about the plotting of the picture.
External procedure DAXIS(integer id; integer array frame;
string xlabel,ylabel; real array params);
This routine will draw a set of axes for you. It will automatically
compute scaling. You give it the array PARAMS (described below) which
contains the maximum and minimum values of X and Y for the function
to be displayed, and DAXIS computes the (linear) transformation from
function coordinates to screen coordinates for you. It also draws
axes, with possibly the XLABEL and YLABEL strings close by if they
are non-null, and puts tick-marks along the axes, and numerical
labels as well as possble (with certain exceptions).
Likewise, the FRAME array contains information about where you want
to plot the picture. This is specified as follows: the screen is
divided, vertically and horizontally, into slices. For instance, if
you ask for 1 slice vertically and 1 slice horizontally, you get the
entire screen. If you ask for 2 slices vertically and 1 slice horizontally,
you will get two rectangles such that the width of each will be the
full screen width and the height of each will be one-half the screen
height. These slices are numbered from 1, starting at the minimum
coordinate. For instance, if you specify 4 vertical slices, the bottom-most
slice will be numbered 1, the next 2, the next 3, and the top-most 4.
In the FRAME array, you give 4 numbers, which are the number of slices
in the X and Y directions, and the slice number in each direction that
the current set of axes is to go in.
You can also specify the rectangle on the screen in terms of absolute
screen coordinates (-512 to +511) in both the X and Y directions.
The routine DAXIS only looks at FRAME[1:6] and PARAMS[1:4], and
produces as output PARAMS[5:17]. It sets up many things that are
not used directly in DAXIS but are instead used in other display
routines, such as DPY. The exact meanings of the words of FRAME
and PARAMS are as follows, although we suggest that the casual
reader skip to the text immediately thereafter:
FRAME[1]=NYSLICES, number of slices in vertical direction
or YMAX, maximum in Y direction
FRAME[2]=YSLICEN, slice number of this display
or YMIN, minimum in Y direction
FRAME[3]=NXSLICES, number of X slices
or XMAX, maximum in X direction
FRAME[4]=XSLICEN, slice number of this display
or XMIN, minimum in X direction
FRAME[5]=TYPE = 0 ⊃ Sample the function if overflow
TYPE = 1 ⊃ Average the function
TYPE = 2 ⊃ Take Envelope of function, FRAME[8] has # pts
TYPE = 3 ⊃ Do direct display
bits 30-32 specify style
30-32 = 0 (0) ⊃ Connect points with linear segment
30-32 = 1 ('10) ⊃ Make step function
30-32 = 2 ('20) ⊃ Just dots
bit 29 neq 0 ⊃ Take absolute value of signal
bit 28 is set to zero by DAXIS, implies `first point' to
DPY. For overlaying, must set this to zero each overlay.
bit 27 neq 0 ⊃ use (xmax,xmin) rather than (xnslices,xslicen)
bit 26 neq 0 ⊃ use (ymax,ymin) rather than (ynslices,yslicen)
bit 25 neq 0 ⊃ inhibit plotting X axis
bit 24 neq 0 ⊃ inhibit plotting Y axis
FRAME[6]=NPTS, for enveloping, gives max window
FRAME[7]=N, Total number of points to be displayed
PARAMS[1]=VXMIN Limits on the Values of the functions
PARAMS[2]=VXMAX
PARAMS[3]=VYMIN
PARAMS[4]=VYMAX
PARAMS[5]=XOFFSET (returned by DAXIS)
PARAMS[6]=XSCALE "
PARAMS[7]=YOFFSET "
PARAMS[8]=YSCALE "
PARAMS[9]=XPOS " (set to VXMIN initially)
PARAMS[10]=TYPE " (copied from FRAME[5])
PARAMS[11]=CSUM " (current sum for averaging mode)
PARAMS[12]=CWIDTH " (number of points left this window)
PARAMS[13]=CDEL " (Increment in value for each X-sample)
PARAMS[14]=NPSUM " (copied from FRAME[6])
PARAMS[15]=XLC " (set to VXMIN initially)
PARAMS[16]=LASTX " (Set to last X pos in screen coords)
PARAMS[17]=LASTY " (Set to last Y pos in screen coords)
Most of this material is included for your information if you ever have
to go inside these arrays. Fortunately, for the more casual user, there
are display routines that do this automatically. This is typified by
the following routine:
External procedure DPYSL(integer id; real array X;
integer N,Nslices,sliceN; string xlabel,ylabel;
real vxmin,vxmax);
This routine takes a real array (X) of N points and displays it. It divides
the screen into NSLICES horizontal strips and places our picture in strip
number SLICEN. For example (to repeat what we mentioned above), setting
NSLICES=SLICEN=1 will give a full-screen image. NSLICES=2 will divide
the screen into upper and lower halves such that SLICEN=1 will be on bottom
and SLICEN=2 will be on top. The X axis will be labeled with the string
XLABEL and the Y axis (vertical axis) will be labeled with the string YLABEL.
The routine will number the vertical axis corresponding to the values in
the X array. It will find, all by itself, the maximum and minimum values
and will display the entire range of the X array in the desired slice.
It can not, however, find the maximum and minimum values of the horizontal
axis, because it only has the index number between 1 and N. For this reason,
we include two more numbers in the calling sequence, which are the maximum
and minimum values to be displayed on the horizontal axis. The number VXMIN
will correspond to the left side (index = 1) and the other number VXMAX
will correspond to the right side (index = N).
An example of the usage of this might be as follows:
! This is extracted from the program JAMA:POKE.SAI[MA,JAM] ;
preload_with -1.19685,0.19685;
real array D[1:2]; ! Windowing function ;
! We assume here that there is a sound file ;
while true do ! already opened on channel 1 in mode '17 ;
begin "HAIRY"
inreal(time,"Starting time"); ! Ask user for time in sound file;
inreal(window,"Window width"); ! And width in seconds;
bgss←time*clock; ! Convert time to sample number;
npnts←window*clock;
N2←npnts; ! We will pad by factor of two;
begin "ALAR"
real array X[0:2*npnts-1]; ! Samples go here;
real array Mag,Rl,Im[0:N2]; ! Transform goes here;
integer id; ! Picture id ;
sndin(1,bgss,npnts,H,X,iflag); ! X will already be cleared;
dsetup(1000+3*npnts,id);
dpysl(id,X,npnts,3,3,time,time+window,"TIME","AMP");
! 3 horizontal strips, signal at top;
weight(X,npnts,D,2); ! Take FFT of window;
artran(X,Rl,Im,N2,false);
for i←0 step 1 until N2 do
Mag[i]←sqrt(Rl[i]↑2+Im[i]↑2); ! Magnitude ;
dpysl(id,Mag,N2,3,2,0,clock/2,"FREQUENCY","MAG");
write(id,0); ! Display spectrum in middle;
drels(id);
end "ALAR";
end "HAIRY";
This program will loop indefinitely, each time asking the user where in
the sound file he wants to display things. It will display the signal at
the top of the sound file and the spectrum below. As often as not, we
prefer to look at the log-magnitude spectrum. this can be done by the
following change to the above program:
artran(X,Rl,Im,N2,false); ! Just like before;
for i←0 step 1 until N2 do
Mag[i]←sqrt(Rl[i]↑2+Im[i]↑2); ! Magnitude ;
xmax←0;
for i←0 step 1 until N2 do
xmax←xmax max Mag[i]; ! Get maximum of spectrum;
for i←0 step 1 until N2 do
if Mag[i] ≤ 0.001*xmax then Mag[i] ← -60.0
else Mag[i] ← 20*log(Mag[i]/xmax)/log(10);
! Convert to decibels, limited to 0 to -60 dB ;
dpysl(id,Mag,N2,3,2,0,clock/2,"FREQUENCY","LOG MAG IN DB");
As soon as people start using the display package, they immediately want
to start doing overlays, or adding things to the images, such as little
pointers or squares or other such things. At first thought, you might
imagine that you can do overlays just by calling DPYSL twice with the
same slice number. Indeed, you will get two displays, one on top of the
other, but they will both have a set of axes, and the axes will not
necessarily line up because DPYSL scales each one vertically depending
on the maximum and minimum values found in the X array. If these maximum
and minimum values are different in the two overlays (and presumably they
will be), then the scaling, and thus the labeling of the axes, will be
different. This makes a bit of a mess on the screen. The correct way
to do it is to use the PARAMS and FRAME array that DPYSL has already
set up with the correct scaling. As an example, we repeat the inner
loop of the previous program, but this time, we will overlay the spectrum
with the unwindowed spectrum, to see the difference:
External real params; ! These is the array that DPYSL uses;
sndin(1,bgss,npnts,H,X,iflag); ! X will already be cleared;
dsetup(1000+3*npnts,id);
dpysl(id,X,npnts,3,3,time,time+window,"TIME","AMP");
! 3 horizontal strips, signal at top;
artran(X,Rl,Im,N2,false); ! Take FFT without window;
for i←0 step 1 until N2 do
Mag[i]←sqrt(Rl[i]↑2+Im[i]↑2); ! Magnitude ;
dpysl(id,Mag,N2,3,2,0,clock/2,"FREQUENCY","MAG");
weight(X,npnts,D,2); ! Now for the windowed FFT;
artran(X,Rl,Im,N2,false);
for i←0 step 1 until N2 do
Mag[i]←sqrt(Rl[i]↑2+Im[i]↑2); ! Magnitude ;
start_code
movei 1,'200; ! Re-initialize parameter array;
movei 2,params; ! (could be done in SAIL with MEMORY construct);
andcam 1,9(2);
move 1,(2); ! Set us back to original X coordinate;
movem 1,8(2);
movem 1,14(2);
end;
dpy(id,params,N2,Mag); ! Create overlay using previous scaling;
write(id,0); ! Display spectrum in middle;
drels(id);
There are at least two pieces of magic here that need explaining
(1) the fact that PARAMS is declared as a REAL rather than a REAL
ARRAY, and (2) the bit with the START_CODE block. They are related.
First off, about the array, the array isn't really an array. Only
SAIL makes arrays, and the display routines are hand-coded in PDP-10
assembly language. Consequently, you cannot access the built-in
FRAME and PARAMS arrays as an array directly in SAIL. There is, however,
a kludge that can be used to get at these data. The first part of
the kludge is to declare the units FRAME and PARAMS as EXTERNAL INTEGER
and EXTERNAL REAL.
The bit with the START_CODE is that the DPY routine is capable of annexing
points to the current graph. That is, instead of calling DPY with all
N2 points at a time, you can give it N2/2 points first, then N2/2 afterwards
in two (or more) separate calls. To accomplish this, DPY keeps certain
key "state" information in locations PARAMS[10] and PARAMS[15]. If you
want your new call on DPY to start at the origin of the display, at the
left of the screen, you will have to reset this state information. Otherwise,
the image will be concatenated onto the end of the current curve. This is
done by resetting the '200 bit in the type word, and by resetting the previous
X coordinate to the origin coordinate. The START_CODE section does exactly
that.
The only problem here is that the scaling will be the scaling of the first
graph displayed. To make this work right, you should display the graph with
the widest dynamic range first. Otherwise, you will get points off the
screen.
The next thing that people typically want to do is to point to things
on the screen, or put boxes around them. This you can do by using the
scaling already present in the PARAMS array. To be exact, you can
define the following things:
real xoffset,xscale,yoffset,yscale,vymin,vymax;
external real params;
! This transforms a value from VXMIN to VXMIN to screen
coordinates (-512 to 511) ;
real procedure TX (real xval);
return(xval * xscale + xoffset);
! This transforms a value from VYMIN to VYMIN to screen
coordinates (-512 to 511) ;
real procedure TY (real yval);
return(yval * yscale + yoffset);
! . . . ommitted part of program that does DPYSL here . . . ;
! Now set up variables for TX and TY to use ;
xoffset ← memory [location(params) + 4, real];
xscale ← memory [location(params) + 5, real];
yoffset ← memory [location(params) + 6, real];
yscale ← memory [location(params) + 7, real];
vymax ← memory [location(params) + 3, real];
vymin ← memory [location(params) + 2, real];
At this point, you can now do display using AIVECT and AVECT and friends
to position the beam using the coordinates of the data. For example, to
position something to X1, Y1, where X1 is in terms of the numbers along
the horizontal axis and Y1 is in terms of the numbers along the vertical
axis (this means VXMIN≤X1≤VXMAX and VYMIN≤Y1≤VYMAX), you could say
AIVECT(id,tx(X1),ty(Y1));
The TX and TY subroutines will transform the data coordinates to screen
coordinates. Obviously, you can freely intermix screen and data coordinates
like this. Thus for instance, to put a little arrow (←) next to some point,
(X1,Y1), we might use the following sequence:
AIVECT(id,tx(X1)+4,ty(Y1)-6); ! To the right and down a little;
DTEXT(id,"←",0,0);
Likewise, we could get a vertical arrow by asking for a 90-degree rotation
of the down-arrow.
The following routine, PNTPLT, is for making fancy multi-overlayed graphs
with circles and dots and stuff like that. The routine was written by
John Gordon of Stanford, and since I don't understand it myself, you will
have to be content with his own explaination, which follows:
External procedure PNTPLT(integer id; boolean new,polar; real array parms,x,y;
integer array fraam,chars;
integer npnts;boolean single;
integer chr,connect;
reference real scale;
string xlabel,ylabel;
boolean quad);
You use this routine in place of DPYSL - that is, you claim a buffer
with DSETUP, and you display the picture afterwards with WRITE, then
release the buffer with DRELS.
If NEW is true, new axes are computed and drawn according to information
in the PARMS and FRAAM arrays. NOTE: THIS ROUTINE MAKES ITS OWN DAXIS
CALL, AND HENCE SETS UP AND READS THE FRAME AND PARAMS ARRAYS THAT ARE
INTERNAL TO THIS PACKAGE. BUT TO LET THE USER SPECIFY HIS OWN NUMBERS IN
THESE ARRAYS, IT MAKES THE USE OF ARRAYS AS PARAMETERS.
PARMS needs to be a one-dimensional, 8-element array. You should load it,
in order, with these 4 real numbers: XMIN, XMAX, YMIN, YMAX. These are
the lower and upper bounds on the values that your point arrays will
assume. PARMS[5:8] are used by this procedure, and you shouldn't touch
them. They're especially important if you want to overlay a new plot on
top of an old one.
FRAAM is a one-dimensional, 5-element array. The first 4 represent screen
coordinate values, and are interpreted according to the bits set in the
5th word of the array. Setting bits mean the following:
'0400 FRAAM[3] & FRAAM[4] are interpreted as actual screen
coordinates, rather than slice numbers (see below).
'1000 Same as '400, except for FRAAM[1] & FRAAM[2]. (Y values)
'2000 Inhibit plotting of the X axis
'4000 Inhibit plotting of the Y axis
All other bits of the word are ignored. The number of slices refers to
how many rectangles you want to cut the screen up into. The number is not
how many cuts are made, but rather, how many pieces you get. If bit '0400
is off, FRAAM[3] should contain the number of slices in the X direction,
and FRAAM[4] should have the actual slice number desired (slices are
numbered from left to right). Likewise for FRAAM[1] and FRAAM[2] for the
Y direction, slices being numbered from bottom to top. If you wish to set
your own rectangle size and position, set bit '400 or '1000 or both in
FRAAM[5], and load FRAAM[1:4] with: Max Y coord, Min Y coord, Max X
coord, Min X coord. Screen coordinates run from -512 to +511 in both
directions.
X is an array of X coordinates, and Y contains the corresponding Y
coordinates, and they should all be within the corresponding bounds given
in the PARMS array.
If POLAR is true, then the X array has magnitudes, and the Y array has
angles in radians, for polar coordinates. You have to set up the PARMS
array with appropriate values (probably + & - the largest magnitude for
both horizontal and vertical.
There are two different ways to specify what characters you wish to use to
represent the points you want to plot. If the same character is used for
all points, then make SINGLE true, and set CHR to the character value
(remember, CHR is an integer, not a string. This means you can ask for
wierd characters like formfeed and altmode, and also default characters).
If CHR is ≥128 (decimal), then default characters will be used:
'200 = 128 ≡ empty square
'201 = 129 ≡ empty triangle
'202 = 130 ≡ empty circle
'203 = 131 ≡ filled-in square
'204 = 132 ≡ filled-in triangle
'205 = 133 ≡ filled-in circle
If you want different characters for each point, you can put these into
CHARS and set SINGLE to false.
NPNTS is the number of points in the X, Y, and CHARS arrays.
CONNECT's value determines how to connect the points:
0 ≡ no connecting lines
1 ≡ connect with solid lines
2 ≡ dotted lines
3 ≡ dased lines
4 ≡ dotted-dashes lines
If you don't like the way this routine computes the scaling factor for
text, you can set your own scale. SCALE=0 means you want automatic scale
computing. SCALE<0 means you want the default size (which is the size of
characters used for everyday use). SCALE>0 means SCALE contains the actual
scaling size you desire.
This routine will plot the X=0 and Y=0 lines for you if you want (i.e,
divide your plot into quadrants), and does so if QUAD is true.
Finally, if you want to label your axes, then put the strings in XLABEL
and YLABEL. If you use slice numbers, the routines allocate space for the
axes and labels if you have any. If you're using your own coordinate
values, then you may be interested in knowing that 35 points are reserved
for each axis, 45 for X labels, and 100 for Y labels. Thus, if you want
to slam your picture against the left margin, you ought to use -477 as the
Min X.
DISPLAY BUFFER FORMAT
This part is for more advanced usage. The novice can skip over the rest
of this section.
The ID word is a SAIL integer which is the address of the first word of
of the buffer that is built by DSETUP. The first word is precisely an
AOBJN pointer to the next free word of the buffer. To be more precise,
the left half of the first word (the word pointed to by ID) is the negative
of the number of remaining words in the buffer, and the right half of the
first word is the absolute address of the last word written into the buffer.
The following words are the buffer itself. To be more explicit:
ID → -N ,, <adr of last word>
<dpy word 1>
<dpy word 2>
. . .
<dpy word M>
< N free words >
where ",," is taken to mean <left halfword>,,<right halfword>.
We can access this through SAIL by use of the MEMORY and LOCATION operations,
like for instance MEMORY[ID,INTEGER] gets you the pointer word (-N,,adr).
MEMORY[ID+1,INTEGER] gets you the first display word, and so on.
When the buffer is written out on the disk by the DWRITE subroutine, the
format is a bit different. It is as follows:
WD CONTENTS COMMENT
---------------------------------------------------------------
0 <ignore> Historical artifact - please ignore
1 Nwds-1 Number of display words minus 1
2 <dpy word 1> First display word
. . .
Nwds+1 <dpy word NWDS> Last display word
Thus, you can read a plot file into core and display it as follows:
! Open PLT file here in mode '10, buffered binary mode
I assume it is open on channel CHAN;
wordin(Chan); ! Flush first word;
Npnts←wordin(Chan)+1; ! This is total word count;
begin "ALAR"
integer array buf[0:Npnts];
arryin(Chan,buf[1],Npnts); ! Read in display;
buf[0]←location(buf[Npnts]);
id←memory[location(buf),integer]; ! Show zero remaining words;
write(id,0); ! Post the display;
end "ALAR"; ! Do not use DRELS - exiting block will delete array;
If you feel you have to be able to use DRELS, then you can do the following:
! Open PLT file here in mode '10, buffered binary mode.
I assume that it is open on channel CHAN ;
wordin(Chan); ! Flush first word;
Npnts←wordin(Chan)+1; ! This is total word count;
dsetup(Npnts,id); ! Allocate a buffer of the proper size;
arryin(Chan,memory[id+1],Npnts);
! Read into data portion of buffer;
memory[id,integer]←id+Npnts;
! Address of last word in buffer, 0 left half;
Now you may use DRELS to release the buffer, since it was claimed by
DSETUP just like normal.
If you are reading in this picture, but you wish to use a larger buffer because
you are going to add some more display to it, you might wish to use the
following code fragment:
! Open PLT file here in mode '10, buffered binary mode.
I assume that it is open on channel CHAN ;
wordin(Chan); ! Flush first word;
Npnts←wordin(Chan)+1; ! This is total word count;
dsetup(Npnts+Np,id); ! Allocate a buffer with NP extra words;
arryin(Chan,memory[id+1],Npnts);
! Read into data portion of buffer;
memory[id,integer]←id+Npnts;
! Address of last word in buffer, 0 left half;
dpb(-Np,point(18,memory[id],17));
! Deposit number of remaining free words in left half;
Other things, such as reading in two pictures and concatenating them we will
leave as exercises for the reader.
The only remaining arcana is what the exact format of the buffer is. This
is quite simple. There are three kinds of words: text words, long vectors,
and short vectors. Everything else should be ignored. The text word is
identified by the fact that the low-order bit is non-zero. In this case,
the word will contain 5 ASCII characters in standard PDP-10 format. A null
character (code 0) should be ignored. For the other instructions, the right
4 bits (bits 32:35, '17) contain an opcode, of which two are defined:
OPCODE=6 LONG VECTOR WORD
BIT OCTAL COMMENT
----------------------------------------------------------------------
0:10 777600,,0 X-value, 2's complement integer
11:21 377,,740000 Y-value
22:24 034000 Brightness. 0 means no change, 2 is normal.
25:27 003400 Character size. 0 means no change, 2 is normal.
29 000100 Mode. 0=relative vector, 1=absolute vector
30:31 000060 Type. 0=visible, 1 and 3=endpoint, 2=invisible
32:35 000017 Opcode = 6
The screen is considered (as we mentioned before) to be 1024 points by
1024 points with the origin (0,0) at the center. The MODE field determines
whether the vector specified is the actual position to move to (absolute
mode) or is the amount to be added into the current (X,Y) coordinates.
If you notice carefully, the coordinates in the long vector word are in
fact 11 bits long, whereas the addressing for the screen is only 10 bits.
This is because the maximum length vector that goes from (-512,-512) to
(511,511) is specified by a difference of (1023,1023), which is an 11-bit
signed quantity. To be able to store differences equal to twice the coordinate
space, there is an extra bit in each coordinate.
The character size refers to any text words that follow. The brightness refers
to all display that follows.
OPCODE=2 SHORT VECTOR WORD
BIT OCTAL COMMENT
----------------------------------------------------------------------
0:6 774000,,0 X-difference, 7-bit 2's complement integer
7:13 003760,,0 Y-difference, 7-bit 2's complement integer
14:15 000014,,0 Type. 0=visible, 1 and 3=endpoint, 2=invisible
16:22 000003,,760000 X-difference
23:29 177000 Y-difference
30:31 000060 Type, as before
A short vector word contains two relative vectors, the coordinates of
which are in the range -64 to +63 (a 7-bit 2's complement number). Each
vector has its own TYPE byte, saying whether it is visible, invisible,
or endpoint.
MISCELANEOUS - INREAL, ININT, DEFAULT
Source: JAMA:IOPACK.SAI[LIB,JAM]
There are three routines that we have already used that perhaps should
be explained at this point:
External procedure INREAL(Reference real X; string S);
External procedure ININT(Reference integer X; string S);
Types out string
S&" = "&cvs(X)&" ← " (in the case of ININT)
or S&" = "&cvg(X)&" ← " (in the case of INREAL)
and then waits for the user to type something. If he does type something,
it is converterted to a real number in the case of INREAL or an integer in
the case of ININT and stored in X and the following is typed:
S&" set to "&cvs(X)&crlf (in the case of ININT)
or S&" set to "&cvg(X)&crlf (in the case of INREAL)
If he does not type something, X is unchanged, and the string
S&" not changed"&crlf
is typed.
External string procedure DEFAULT(string Inp; reference string Dev;
string Devext,Devdev);
File name scanner. File names may be typed in any of the
following forms:
dev:name.ext[pp,pn]
dev:name.ext[pp,pn
dev:name.ext[pp,
dev:name.ext[pp
dev:name.ext[
dev:name.ext
dev:name
name
If the PN is ommitted, the PN of the current alias is used. If the entire
PPN is ommitted, the PPN of the current alias is used. If the extension
is ommitted, DEFEXT will be substituted. If the device is ommitted, DEFDEV
will be used. The procedure returns the filename-extention-ppn string as
its value, ready for a SAIL LOOKUP or ENTER, and returns the device string
in DEV (without colon), all ready for a SAIL OPEN. It puts in the right bracket
and comma for the ppn automatically even if you leave it off. Thus, you
can get and type out the name of the file with the following sequence:
outstr("Input MSB file: ");
sfile←default(inchwl,sdev,"MSB","DSKM");
outstr("Reading "&sdev":"&sfile&crlf);
MISCELANEOUS - LINE SEGMENT FUNCTION IO
Numerous programs use line segment functions to specify various
things. For this purpose, there are routines to read and write
line segment functions in a format compatible for the speech processing
routines as well as for the MUSIC program. The format for a line
segment function is the following:
SEG(<name1>);
<val 1> <time 1>
<val 2> <time 2>
. . .
<val N> <time N>
SEG(<name2>);
<val 1> <time 1>
<val 2> <time 2>
. . .
<val M> <time M>
As a more concrete example, we might give the following:
SEG(AMP1);
0 0
1.0 0.32
1.0 0.75
0 1.20
SEG(AMP2);
1 0
0.5 0.5
0.25 1.0
0.125 1.5
0.0625 2.0
0.03125 2.5
0 3.0
This example shows two functions, one named AMP1 and the other named
AMP2. The first one is 1.20 seconds long and shows a function that
starts at zero, goes to one, stays at one for a while, then comes down.
The second function starts at one and decays roughly exponentially, then
is truncated to zero over a 3 second period. Note that the first column
in each case is the value of the function, and the second column is the
independent variable (time in this case).
Functions in the MUSIC program must be defined with the independent
variable running strictly from 1 to 100.0, so the above two functions
would have to be defined, for feeding to MUSIC, as follows:
SEG(AMP1);
0 1
1.0 27.4
1.0 62.875
0 100.0
SEG(AMP2);
1 1
0.5 17.5
0.25 34
0.125 50.5
0.0625 67
0.03125 83.5
0 100
In any case, the subroutines will accept them either way.
The routine that reads in a list of functions is the following:
External boolean procedure SEGSYN(string dev,file;
reference record_pointer(any_class) L; reference integer nfns);
This routine opens the disk file, using GETCHAN to obtain the IO channel
number. It skips over the E directory page, if there is any, and reads
in all the functions in the file. These are loaded into SAIL records
of the following class:
record_class seg(integer type; record_pointer(seg) next,last;
string name;
integer nsegs;
real mintime,maxtime,minval,maxval;
real array times,values);
SEGSYN then returns to you the pointer to the record containing the
first function in the list in the var)able L. It will return FALSE
normally (upon success) and the number of line segment functions it
read in NFNS. If the file does not exist or is garbled, it returns
TRUE, and L will be NULL_RECORD and NFNS will be set to zero.
The record representing each function will be linked to the next record
in the list and to the previous record in the list by the NEXT and LAST
fields in the record. The integer, NSEGS, gives the number of breakpoints
in the functions. The arrays TIMES and VALUES will be indexed [1:NSEGS].
The real variables MINTIME and MAXTIME are just the minimum and maximum
times in the entire TIME array for this particular function. Likewise,
the MINVAL and MAXVAL variables contain the minimum and maximum values
of the VALUES array. The string variable NAME will be the text name
found for the function (after the keyword "SEG" and the paren in the
original text). This name will be made all majescules, and all spaces
and tabs eliminated. The TYPE variable will allways be zero - it is
reserved for future expansion.
If you wish, for example, to read in several files and concatenate their
function lists to make one big function list (as is done in the program
LPS, for example), the following program fragment may be used to do that:
record_pointer(seg) top,newf,r;
top←null_record;
while true do
begin "FNRD"
string wfile,wdev;
outstr("Warp function file (CR to xit): ");
wfile←default(s←inchwl,wdev,"FUN","DSK");
if length(s)≤0 ∨ s='175 then done "FNRD";
if ¬segsyn(wdev,wfile,newf,i) then
begin "NF"
r←newf; ! Go down to the end of this list;
while seg:next[r] neq null_record do r←seg:next[r];
seg:next[r]←top; ! Concatenate previous list to this one;
if top neq null_record then
seg:last[top]←r;
top←newf; ! Make us the new top of the list;
end "NF"
else outstr("File "&wdev&":"&wfile&" not found"&crlf);
end "FNRD";
The following subroutine will print out the names of all the functions
in a list:
Comment PRINT_LIST - Prints out the function names;
Procedure PRINT_LIST(record_pointer(any_class) r);
begin
integer cc;
string ts;
boolean fst;
fst←true;
cc←0;
while r neq null_record do
begin "PR"
define nsp="5";
ts←seg:name[r]&"(";
if seg:type[r]=0 then ts←ts&"SEG["&cvs(seg:nsegs[r])&"])"
else if seg:type[r]=1 then ts←ts&"SSEG["&cvs(seg:nsegs[r])&"])"
else if seg:type[r]=2 then ts←ts&"SYNTH)";
if cc+length(ts)>80 then
begin "EOL"
outstr(crlf);
cc←length(ts);
end "EOL"
else begin "SL"
cc←cc+nsp+length(ts);
if fst then fst←false
else outstr(" ");
end "SL";
outstr(ts);
r←seg:next[r];
end "PR";
outstr(crlf);
end;
That routine will print out the name of each function, its type (currently
only SEG will be found) and the number of breakpoints present.
Once you have read in your function list, you might like to locate a
certain one. For instance, you might want to ask the user to type the
name of the function he wants, then search for it in the list. The following
code fragment will do that:
outstr("Function name: ");
setbreak(1,"",'40&'11,"INKS"); ! No break, make all capitals;
tmpstr←inchwl;
s←scan(tmpstr,1,brk); ! Make all caps, kill spaces;
r←top;
while r neq null_record do ! Go down list looking for it;
begin "FNDT"
if equ(s,seg:name[r]) then done "FNDT";
r←seg:next[r];
end "FNDT";
if r neq null_record then
begin "FNDTR"
! **** Your function is now in R, do with it as you will **** ;
end "FNDTR"
else outstr("Function "&s&" not in list"&crlf);
Once you have this function in memory and isolated, you will probably want
to read it sequentially, interpolating linearly between adjacent breakpoints.
For a strictly sequential reading, the following pair of subroutines will
do that (though not in the most efficient manner imaginable):
! SETFSCN, FSCN - function scanners;
define tscale="D[1]",toffset="D[2]", Comment - Converts input index into time;
vsc="D[3]",voff="D[4]", Comment - Converts value into output range;
ind="memory[location(D[5]),integer]",
vscale="D[6]",voffset="D[7]"; Comment - Converts time into value;
procedure SETFSCN(record_pointer(seg) fn;
real outlo,outhi,inlo,inhi;
real array D);
begin
tscale←(seg:maxtime[fn]-seg:mintime[fn])/(inhi-inlo);
toffset←-inlo*tscale+seg:mintime[fn];
vsc←(outhi-outlo)/(seg:maxval[fn]-seg:minval[fn]);
voff←outlo-vsc*seg:minval[fn];
ind←0;
end;
Real procedure FSCN(real inp; record_pointer(seg) fn; real array D);
begin
real ltime,vl;
ltime←inp*tscale+toffset;
if ind=0 ∨ (ind<seg:nsegs[fn] ∧
ltime≥seg:times[fn][ind+1])
then begin "SETI"
ind←ind+1;
if ind<seg:nsegs[fn] then
begin "NEWV"
vscale←(seg:values[fn][ind+1]-seg:values[fn][ind])/
(seg:times[fn][ind+1]-seg:times[fn][ind]);
voffset←-vscale*seg:times[fn][ind]+seg:values[fn][ind];
end "NEWV";
end "SETI";
if ind<seg:nsegs[fn] then
vl←ltime*vscale+voffset
else vl←seg:values[fn][ind];
return(vsc*vl+voff);
end;
Since these are among the more opaque routines ever written, a few words of
commentary seem to be in order. The process of reading off a function is
done in two parts: setup is first (SETFSCN). It loads up a little array
(that you must provide - 1 for each function being read at the same time).
The second part is delivering each point of the function sequentially. There
are two scalings involved. The function itself will be defined over some time
and will take on some range of values. These are, as often as not, not exactly
the terms in which you wish to access them. For instance, if the abcissa is
defined in terms of time, you may want it in sample numbers. The values may
be in frequency, but you may want them in multiplied by MAG already. To do this,
you give the ranges for the abcissa and the ordinate separately to the SETFSCN
routine. It will scale both the time and value axes to correspond to your
limits. For example, if you wanted to index the successive values of the
function as if it were an array that went from 1 to 1000, and you wanted the
output values to vary from 0 to 1, inclusive, you might do the following:
real array Darr[1:7];
setfscn(func,0,1.0,1,1000,Darr); ! Load scaling into DARR;
for i←1 step 1 until 1000 do ! Step through the promised range;
begin
real val;
val←fscn(i,func,Darr); ! Read off next value;
! *** Do what you will with this value *** ;
end;
Regardless of how the function was defined, it will then read off sequential
values from 0 to 1.0 over the index range of 1 to 1000. The function could
have been defined any way at all, like either form of AMP1 or AMP2 that we
showed above, and the values delivered will be the same. If you want to preserve
the original range of the function, whatever it may be, the following
setup will do that:
setfscn(func,seg:minval[func],seg:maxval[func],1,1000,Darr);
This causes the numbers read out to be between the minimum and maximum values
that were present in the function definition to be used.
If you want to create your own functions inside your program and link them
up one way or the other, you can do something like the following:
procedure NEWFUN(real array ntimes,nvalues; integer nsegs;
reference record_pointer(seg) top);
begin
External integer procedure ARMAK(integer lb,ub,Ndims);
record_pointer(seg) r;
integer j;
r←new_record(seg);
seg:type[r]←0;
seg:nsegs[r]←nsegs;
memory[location(seg:times[r])] ← armak(1,nsegs,1);
memory[location(seg:values[r])] ← armak(1,nsegs,1);
seg:maxtime[r]←seg:mintime[r]←ntimes[1];
seg:maxval[r]←seg:minval[r]←nvalues[1];
for j←1 step 1 until nsegs do
begin "XFR"
seg:times[r][j]←ntimes[j];
seg:values[r][j]←nvalues[j];
seg:maxtime[r]←seg:maxtime[r] max ntimes[j];
seg:mintime[r]←seg:mintime[r] min ntimes[j];
seg:maxval[r]←seg:maxval[r] max nvalues[j];
seg:minval[r]←seg:minval[r] min nvalues[j];
end "XFR";
seg:next[r]←top;
if top neq null_record then
seg:last[top]←r;
top←r;
end;
This copies the NTIMES and NVALUES arrays into new arrays in a new record,
then places that new record at the top of the list.
External boolean procedure SSWRT(string dev,file;
reference record_pointer(any_class) L);
This procedure will write out a list of functions on the disk as a text
file. There is nothing magical about it, so we will not explain it any
further here. It gets its IO channel number via GETCHAN. It returns
FALSE upon success and TRUE if there was any error (that is, if it got
the error return from the ENTER on the file name you gave).